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1  Abstract 
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The  irony  embedded  in  optic  tow  computation  isthat  despite  its  not  having  been  studied  mathematically 
until  about  20  years  ago,  humans  have  been  performing  the  related  calculations  innately  since  their  creation. 
T  he  process  can  be  thought  of  as  ..lling  in  the  gaps  that  occur  between  the  images  we  see  changing  in  time. 
For  the  human  mind,  the  transition  is  automatic  from  what  an  image  looks  like  at  one  instant  to  what  it 
looks  like  at  the  next.  Our  brain  intuitively  generates  a  continuous  tow  of  what  it  perceives  is  being  viewed 
by  the  eyes. 

Manifested  in  a  mathematical  sense,  optic  tow  serves  as  a  mechanism  to  describe  the  movement  of 
objects  in  a  digital  image  sequence  using  a  tow  ..eld  of  grayscale  intensity.  This  particular  study  has  focused 
upon  developing  a  usable  mathematical  implementation  for  two  of  the  various  algorithms  which  have  been 
proposed  to  compute  optic  tow.  One  method  relies  upon  calculus  of  variations  to  regularize  the  problem 
while  the  other  utilizes  a  wavelet  based  multi-scale  approach. 

Possible  applications  for  such  research  range  from  data  compression  of  video  sequences  to  the  develop¬ 
ment  of  a  moreep  dent  way  to  analyze  the  time  varying  one-dimensional  acoustic  imagery  supplied  by  sonar 
systems. 
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x  =  (a;  1,0:2)  Cartesian  coordinate  of  a  pixel  in  a  2-D  image 
/(x;  t)  G  rayscale  intensity  of  the  image  at  position  x  and  time  t 

d 

— —  Partial  derivative  with  respect  to  xr  where  i  =  1,2 

OXi 

d 

—  Partial  derivative  with  respect  to  t 

at 

v(x; t)  =  (vi,v2)  =  Optic  Flow  vector 

(  dl  dl\ 

V/  =(——  - —  Gradient  of/ 

\axi  0x2  J 

g{x)  Complex  Conjugate  of  the  function  g 


{fig)  =  f(x)g(x)dxidx2  Inner  Product  of  /  and  g 


Generic  symbol  for  a  1-D  wavelet 
<j>  Generic  symbol  for  a  1-D  scaling  function 
9n  2-D  mother  wavelet  with  n  =  1,  ..,7V  where  TV  =  3  or  4 
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=  ip(xi)(j>(x2) 

High-Low  component 

04 

i- I 

II 

Low-High  component 
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=  4>(xi)ip(x2) 

High-High  Component 
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04 
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i- I 

II 

Low-Low  Component 

=  r(x  -  u) 

where  u  is  a  2-D  continuous  translation  index 

CM 

=  si  xen 

—  \  where  s  is  a  continuous  scaling  index 

T  his  forms  a  2-D  continuous  mother  wavelet  family 
9jk (x)  =  2 (2jx  -  k)  where  j  e  Z  and  k  e  Z2 

This  forms  a  2-D  discrete  mother  wavelet  family 
where  ^"k  is  located  around  (2'  jki,2'  jk2)  and  spreads 
over  a  domain  of  size  proportional  to  2'  ■» 

L2{ R)  Set  of  all  square-integrabie  functions  in  R 
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4.1  Introduction 

What  is  it  about  motion  that  human  beings  ..nd  so  interesting?  Why  are  we  so  much  more  likely  to 
follow  a  moving  object  with  our  eyes  than  we  are  to  let  them  rest  on  something  which  is  static?  Perhaps  it 
is  because  of  the  questions  we  learn  to  automatically  ask  ourselves:  Where  is  the  object  going?  W  here  is  it 
comingfrom?  How  fast  is  it  moving?  Do  I  need  to  move  as  well  in  order  to  avoid  its  path?  A nswers  to  this 
cycle  of  questions  must  continually  be  found  as  one  carries  out  the  process  of  evading  collision.  Perhaps  the 
intriguing  factor  about  motion  is  that  it  always  takes  us  to  someplace  new,  someplace  dinerent  from  where 
we  just  were.  T  he  old  Italian  proverb  speaks  truly,  "A  rolling  stone  gathers  no  moss."  As  humans  we  yearn 
to  see  what  we  have  not  yet  seen,  and  only  through  movement  may  we  accomplish  this  goal  and  turn  our 
view  to  the  unexplored. 

This  same  fascination  with  motion  carries  over  to  the  use  of  a  computer  in  our  daily  lives.  Namely,  it 
is  while  watching  video  that  our  minds  are  most  captivated  by  the  screen.  Neither  text,  nor  graphics,  nor 
sound  can  convey  so  much  information  to  the  brain  so  quickly.  Video  is  essentially  the  playing  of  image 
sequences  and  it  is  only  here  that  one  may  watch  objects  moving  through  space  in  time.  But  how  does  the 
computer  view  this  wealth  of  data?  Does  it  assign  any  signi..cance  to  what  is  happening  in  the  video? 

Theanswer,  put  quitesimply,  is  no.  A  computer  devotes  absolutely  no  enort  to  understanding  how  the 
objects  are  moving  in  a  digital  image  sequence.  It  merely  displays  frame  after  frame  and  leaves  it  up  to  the 
user  to  mentally  ..II  in  the  gaps  between  images  in  order  to  perceive  smooth  motion.  For  instance,  let's  say 
a  video  clip  shows  a  baseball  being  thrown  from  one  person  to  another.  E  ven  though  the  ball  is  at  a  certain 
position  in  one  image  frame  and  at  another  position  in  the  next,  the  computer  has  no  way  of  knowing  where 
the  ball  probably  was  for  the  time  in  between  the  two  images.  T  he  human  mind,  on  the  other  hand,  is 
intelligent  enough  to  ..gure  out  that  the  ball  must  have  been  somewhere  between  those  two  positions  during 
that  time.  T  hus  our  brain  creates  a  continuous  fow  of  motion  instead  of  seeing  just  a  series  of  discrete  still 
images  being  played  in  succession. 

An  interesting  problem  here  results.  Can  we  solve  for  this  fow  in  a  digital  sense?  What  kind  of 
information  would  this  provide?  How  can  such  a  question  be  tied  to  applied  mathematics?  The  answers 
lie  in  the  concept  of  optic  tow.  At  the  most  basic  I  e/el,  it  can  be  understood  as  a  mechanism  for  describing 
the  motion  captured  within  an  image  sequence. 

Theoretical  development  for  this  scheme  has  its  foundation  in  the  seminal  paper  by  Horn  and  Schunck 
[14]  and  has  later  been  expounded  upon  by  the  publications  of  Christophe  Bernard  [3],  [4],  [5],  [6],  [7],  In 
both  cases,  the  authors  do  not  provide  the  code  they  use  to  generate  their  results.  Rather,  they  leave  their 
theories  as  virtual  "black  boxes"  which  may  be  employed  to  extract  information  relating  to  the  motion  of 
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objects  appearing  in  the  image  sequence. 

As  often  will  happen  in  the  undertaking  of  a  project  such  as  this,  more  questions  have  surfaced  than 
were  there  at  the  start.  Many  riddles  still  remain  surrounding  optic  tow.  For  instance,  what  is  the  best 
way  to  solve  for  it?  J  ust  how  sound  are  the  physics  which  constitute  the  theory  behind  optic  tow?  Finally, 
how  useful  is  the  information  gained  when  applied  to  a  practical  context?  Despite  such  issues  still  not 
being  resolved,  new  ground  has  certainly  been  broken  in  the  quest  to  present  optic  tow  as  a  ..eld  of  applied 
mathematics  which  is  open  to  experimentation.  A  concrete  path  for  implementing  the  idea  has  ..nally  been 
laid  down. 

4.2  Outline  of  Approach 

The  body  of  this  paper  is  divided  into  ..ve  main  sections.  T  he  ..rst  serves  to  illustrate  the  details  of 
optic  tow,  how  the  concept  is  applied  to  an  image  sequence,  what  kind  of  equation  is  being  solved,  and 
exactly  what  information  is  being  extracted.  As  our  study  has  taken  two  completely  dinerent  directions 
in  trying  to  solve  for  this  optic  tow,  the  next  two  sections  of  the  paper  are  devoted  to  an  explanation  of 
each  method.  At  the  outset  of  the  project,  we  sought  to  treat  the  problem  using  a  wavelet  based  multiscale 
approach  in  the  same  fashion  as  Bernard  [3],  [4],  [5],  [6],  [7],  The  actual  implementation  required  to  do 
this  proved  itself  much  more  complicated  than  originally  envisioned.  G  iving  up  on  the  problem  was  not  an 
option,  however,  and  so  the  decision  was  made  to  begin  a  new  approach  in  parallel  to  our  ongoing  enorts. 
T  his  attempt  relied  upon  the  technique  of  calculus  of  variations  in  a  similar  manner  as  Florn  and  Schunck 
[14],  As  our  research  currently  stands,  the  protocol  for  implementing  each  method  has  been  developed  to  the 
extent  that  they  can  operate  on  simp li . .ed  approximations  to  motion  being  captured  in  an  image  sequence. 
T  he  fourth  section  of  the  paper  evaluates  the  results  attained  using  these  two  approaches.  Conclusions  are 
presented  in  the  ..fth  and  ..nal  section  of  the  paper.  Various  issues  will  be  discussed  here  including:  which 
method  appears  to  work  better;  what  opportunities  exist  to  adjust  each  model  to  better  suit  the  problem; 
what  issues  are  left  unresolved;  what  are  some  of  the  possible  future  applications  for  optic  fow;  and  how 
does  the  overall  idea  behind  optic  fow  stand  against  other  problems  examined  in  applied  mathematics. 

4.3  Desired  Result 

The  theory  behind  this  project  is  both  powerful  and  highly  desirable  for  application  to  a  context  where 
determination  of  an  object's  motion  is  the  foremost  goal.  Without  a  tangible  way  of  utilizing  it,  however, 
a  theory  doesn't  amount  to  having  much  value  at  all.  This  is  especially  true  for  the  Navy  and  the  military 
as  a  whole  which  places  an  indisputable  requirement  on  practicality.  Therefore,  the  primary  focus  of  the 
research  conducted  for  this  project  has  been  to  produce  a  deliverable  which  harnesses  the  potential  of  this 
theory  and  sets  it  down  in  a  usable  code. 
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Asa  second  bene.t  accompanying  our  study,  we  have  performed  a  mathematical  comparative  analysis 
on  thedinerent  methods  for  implementation.  In  each  case,  observations  were  made  regarding  the  advantages 
and  shortcomings  of  that  particular  approach  .  E  nectively,  this  process  may  be  considered  the  optimization 
of  the  tools  we  have  created  for  practical  use.  T  his  is  the  very  crux  of  the  original  project  proposed  because 
it  answers  the  questions  of  how  one  might  take  an  existing  mathematical  theory  and  apply  it.  T  he  knowledge 
gained  throughout  the  development  of  such  implementation  has  been  an  invaluable  reward  to  the  writer  and 
hopefully  it  will  serve  as  the  same  for  those  who  seek  to  build  upon  its  foundations. 


5  Optic  Flow 

5.1  Concept  Introduced 

The  irony  embedded  in  optic  tow  computation  isthat  despite  its  not  having  been  studied  mathematically 
until  close  to  20  years  ago,  human  beings  have  been  performing  the  related  calculations  innately  since  their 
creation.  In  its  most  common  sense,  optic  tow  measurement  is  the  description  of  images  as  they  change  in 
space  and  time.  The  process  can  be  thought  of  as  ..[ling  in  the  gaps  that  occur  between  the  instantaneous 
images  we  see.  For  the  human  mind,  thetransition  is  automatic  from  what  an  image  looks  like  at  one  instant 
to  what  it  looks  like  at  the  next.  Our  brain  intuitively  generates  a  continuous  tow  of  what  it  perceives  is 
being  viewed  by  the  eyes  [4], 

The  digital  analogy  is  not  so  simple.  Image  sequences,  or  video,  is  treated  by  a  computer  as  a  set  of 
discrete  pictures  which  when  cycled  through  quickly,  create  for  the  viewer  the  enect  of  a  scene  changing 
as  time  progresses.  If  the  images  were  not  recorded  at  a  rapid  enough  rate,  when  played  back,  the  video 
will  appear  blurred  and  choppy.  T  his  phenomenon  occurs  because  the  computer  does  not  have  the  inherent 
intelligence  to  say  that  what  was  at  one  place  at  one  time  must  have  gotten  to  its  next  position  along  a 
relatively  smooth  and  logical  path.  It  cannot  automatically  assume  that  for  those  instants  when  it  does  not 
have  a  displayable  image,  the  objects  in  motion  were  at  some  point  in  between  the  recorded  positions.  T  his 
is  exactly  the  problem  which  spawned  the  concept  of  optic  X ow  [4], 

In  one  of  the  most  de..nitive  works  analyzing  the  subject,  C  hristophe  Bernard  describes  optic  fow  as  "a 
visual  displacement  fow  . . eld  that  can  be  used  to  explain  changes  in  an  image  sequence"  [3],  [4],  In  other 
words  it  can  be  thought  of  as  a  vector  ..eld  with  its  arrows  pointing  in  the  direction  of  movement  for  the 
objects  in  motion  within  the  video.  He  provides  a  very  clear  illustration  of  this  in  Figures  1  and  2: 
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Figure  1:  Car  Scene 


Figure  2:  Measured  Optic  Flow 


Flerethree  cars  are  in  motion:  a  light  colored  taxi  making  a  turn  onto  a  side  street,  as  well  as  a  dark  car 
and  a  dark  van  which  are  driving  in  opposing  directions  on  the  main  street.  The  corresponding  measured 
image  tow  graph  shows  a  high  concentration  of  vector  arrows  in  the  location  of  each  of  the  three  moving 
objects  representing  their  direction  of  travel.  In  those  areas  of  the  image  where  there  is  no  motion,  e.g. 
the  parked  cars  and  building,  very  few  arrows  appear  and  none  have  signi . . cant  magnitude.  Our  goal  is  to 
be  able  to  solve  for  these  vectors  based  on  a  mathematical  model  for  what  motion  is  going  on  in  the  image 
sequence. 

Observing  this  vector  ..eld,  we  see  that  the  velocity  of  the  grayscale  intensity  is  being  represented  from 
an  Eulerian  perspective  since  it  gives  us  a  snapshot  of  the  optic  X ow  as  a  function  of  pixel  position.  T  his 
allows  us  to  visualize  how  the  image  sequence  is  changing  at  each  instant  in  time  as  opposed  to  a  Lagrangian 
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approach  which  would  display  a  vector  ..eld  that  would  trace  the  path  of  individual  objects  as  they  moved 
throughout  time  [17,  pg.  304], 

The  immediate  question  that  must  be  asked  is  exactly  how  a  computer  is  supposed  to  correctly  detect  the 
movement  of  objects  in  an  image  sequence.  Surprisingly,  the  answer  is  not  overly  complicated.  Every  image 
stored  digitally  can  be  thought  of  as  a  three-dimensional  surface  with  the  x  and  y  coordinates  describing 
the  location  of  a  pixel  within  the  image  and  the  z  coordinate  relating  to  the  grayscale,  or  intensity,  of  that 
speci..c  pixel.  This  scheme  is  easily  visualized  by  examining  Figures  3  and  4: 


F i gure  4:  Digital  Representation  of  Grayscale  I  mage 
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Within  the  picture  of  the  stapler  sitting  on  the  desk,  there  exist  varying  grayscale  intensities.  These 
directly  correspond  to  thedinerent  heights  of  the  surface  portrayed  in  F  igure  4.  W  here  the  image  in  Figure 
3  is  white,  the  surface  in  Figure  4  is  very  high.  Conversely,  where  the  image  is  black,  the  surface  is  low. 
NOTE:  T  he  color  appearing  in  F  igure  4  serves  merely  to  vi sually  accentuate  the  shape  of  the  surface  and 
has  no  correlation  to  color  appearing  or  not  appeari  ng  in  the  image.  Also,  the  labels  A,  B,  C,  and  D  appear 
only  as  a  guide  for  the  reader  to  ..nd  the  corresponding  corners  between  the  image  and  the  surface. 

As  time  progresses  and  objects  move  within  a  scene,  the  intensities  of  the  individual  pixels  will  invariably 
change  This  movement  is  visualized  by  the  computer  as  the  three-dimensional  surface  changing  shape  It 
is  here  that  mathematics  enters  the  problem  as  surfaces  and  their  changing  shape  in  space  and  time  are 
describable  using  partial  dinerential  equations  (P D  Es). 

5.2  P  rocess  I  llustrated 

The  task  now  becomes  to  ..t  a  mathematical  model  to  the  concept  of  optic  fow.  Our  ..rst  step  in  this 
process  is  to  decide  on  a  set  of  assumptions  which  make  the  problem  more  approachable  at  ..rst  glance.  T  wo 
critical  assumptions  made  by  H orn  &  Schunck  and  Bernard  are:  the  objects  appearing  in  the  image  are 
receiving  uniform  illumination  (i.e.  the  ambient  lighting  is  not  changing);  and  the  grayscale  intensity  of  a 
particular  point  on  the  object  is  constant  over  local  space  and  time  [14],  [3],  [4],  [5],  [6],  [7],  Although  this 
may  seem  like  a  gross  oversimpli.. cation  of  real-life  conditions,  these  constraining  assumptions  actually  make 
sense  when  taken  in  terms  of  the  practicality  of  optic  fow. 

At  its  root,  optic  fow  is  a  mechanism  for  representing  the  motion  of  objects  in  an  image  sequence  using 
vector  arrows.  I  n  order  to  ..nd  these  arrows,  we  choose  a  starting  or  "base"  image  in  the  sequence  in  which 
all  the  objects  pictured  are  said  to  be  in  their  original  position.  As  time  progresses,  the  objects  move  and 
so  their  pixel  positions  change  from  image  frame  to  image  frame.  The  concept  of  optic  fow  is  to  track 
these  moving  objects  and  assign  to  this  motion  a  set  of  vector  arrows  from  frame  to  frame.  As  is  the  case 
with  any  mathematical  model,  there  will  be  some  error  in  optic  fow's  ability  to  predict.  There  will  be  a 
dinerence  between  the  actual  image  frames  and  what  is  predicted  from  the  base  image  and  the  fow  ..elds. 
T  his  dinerence  is  known  as  the  correction  factor.  Since  every  frame  is  known,  the  correction  factor  can  be 
solved  for  exactly  and  so  perfect  reconstruction  can  occur  for  subsequent  images  using  the  base  image,  the 
optic  fow  ..eld,  and  the  correction  factor.  In  essence,  this  says  that  we  can  project  into  the  future  of  the 
image  sequence  from  a  certain  starting  point  in  time. 

Now,  we  ask  what  happens  when  one  of  our  constraining  assumptions  is  violated?  This  could  very 
realistically  happen  when  the  illumination  in  the  image  changes  from  frame  to  frame  and  one  of  the  objects 
in  motion  now  appears  dinerently  in  its  shading/  grayscale  intensity.  We  may  no  longer  use  the  original  base 
image  to  project  forward  in  time  with  optic  fow.  T  he  overall  process,  however,  doesn't  come  to  a  permanent 
halt.  We  merely  begin  over  again  by  setting  the  next  frame  in  the  sequence  to  be  the  new  base  image. 
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From  there  we  begin  tracking  the  motion  of  objects  from  their  new  "original"  position.  F  igure  5  illustrates 
this  process. 


IMAGE  SEQUENCE 


Correction  Correction 

Factor  Factor 
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Field 
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Field 


image  Later  On 
in  Sequence 


Change  in 
illumination 


CHOOSE  NEW  BASE  IMAGE 
AND  REPEAT  PROCESS 


F  i  gu  re  5:  Schemati  c  for  0  pti  c  F I  ow  P  rocess 


There  are  two  ways  of  representing  the  motion  of  objects  smoothly  in  an  image  sequence;  either  use 
more  image  frames  per  unit  time  or  actually  track  the  object  movement  using  vector  arrows  in  optic  tow. 
According  to  the  mdthod  just  described,  optic  tow  provides  the  distinct  advantage  of  using  much  less  digital 
storage  space.  Byte  requirements  for  vector  ..elds /  correction  factors  as  opposed  to  full  images  are  markedly 
less.  Thus  consider  an  average  video  sequence  which  uses  32  frames  per  second.  Even  if  it  is  undergoing 
rapid  illumination  changes  and  a  new  base  image  is  having  to  be  chosen  about  once  a  second,  optic  tow  is 
still  able  to  replace  about  31  images  with  vector  ..elds  [5],  [6], 

5.3  Governing  Equation 

H avin g  stepped  through  that  process,  let  us  now  turn  to  the  development  of  the  governing  partial 
dinerential  equation  (P  DE )  for  optic  tow  [14],  Starting  from  the  constraining  assumption  that  the  grayscale 
intensity  of  an  object  is  the  same  over  local  space  and  time  (and  using  the  notation  previously  de..ned),  we 
arrive  at  the  following  equivalence  relation: 
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I(xi,X2,t)  =  J(x  1  +  Sx  1,X2  +  5x2, t  +  5t). 


(1) 


In  ©<panding  the  right  hand  side  of  Equation  (1)  around  the  point  (xi ,x2,t)  we  see  that 


dl  dl  dl 

I(xi,X2,t)  =  I(xi,x2,t)  +  (5xi^-  +5x  2~^  +St~t  +£- 


(2) 


where  e  is  the  error  term.  If  we  consider  the  limit  of  going  over  the  most  localized  space  and  time,  then 
e  ->■  0  .  The  next  step  is  to  subtract  I(xi,x2,t)  from  both  sides  of  Equation  (2)  and  then  divide  both  sides 
by  5t  .  This  yields 

n  _  Sx-\  dl  |  Sx?  dl  ,  dl 

St  dx\  St  dx2  dt ' 


which  is  equivalent  in  the  prede..ned  notation  to  be 


0 


dl  dl  dl 
~d^Vl+~d^V2  +  -dt 


(3) 


As  Equation  (3)  is  the  governing  P DE  for  optic  tow,  this  represents  the  beginning  of  our  analysis. 

The  image  intensity  function  I  is  known  at  each  pixel  x  at  each  time  t.  Therefore  we  may  state  that 

dl  dl  dl 

the  partial  derivatives  of  the  intensity - , - ,  and  —  can  be  found.  Thus  in  this  governing  equation 

dx i  <9x2  dt 

for  optic  tow,  also  known  as  the  intensity  constraint  equation,  we  see  that  there  are  two  unknowns,  namely 
the  components  of  the  velocity  v\  and  v2.  In  a  single  linear  equation  such  as  this  where  there  are  two 
unknowns,  the  problem  is  said  to  be  ill-posed.  We  cannot  solve  for  both  components  of  the  velocity  without 
an  additional  constraint.  F  igure  6  illustrates  this  concept. 
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Figure  6:  1 1 1- Posed  ness  of  Optic  Flow  Problem 


In  this  example  ..gure,  a  level  curve  of  the  intensity  function  can  be  thought  of  as  a  part  of  the  image 
which  has  the  same  grayscale  intensity  throughout.  By  de. . ni ti on  of  gradient,  the  vector  VI  = 
is  perpendicular  to  the  level  curves.  The  reader  will  see  that  using  the  governing  equation,  one  can  always 
solve  for  the  component  of  the  optic  fow  which  is  in  the  direction  of  the  gradient  of  the  intensity  function 
for  any  given  point  x. 

Equation  (3)  can  be  rewritten, 


/  d I  dl  \ 
\9xi’  dx2  ) 


■  {v\,V2)  =  V /  •  V 


di 

df 


(4) 


If  we  apply  the  d  e.  .ni  ti  on  of  dot  product  that  a  •  b  =  \\a\\  H&ll  cos  <9  where  6  is  the  angle  between  the  two 
vectors,  we  arrive  at  Equation  (5), 


VI- v  =  ||V/||  ||v||cos0 


m 

dt 


(5) 
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Dividing  through  by  ||V/||  gives  us  the  component  of  the  optic  fow  vector  in  the  direction  of  the  velocity 
gradient, 


|| v  ||  cos  6 


dl_ 

_ M. _ 

l|Vi||  ' 


(6) 


This  component  of  the  optic  X ow  vector,  which  is  labeled  in  F  igure  6,  can  always  be  solved  for  using  the 
governing  P  DE .  W  hat  must  be  noted  from  Figure  6,  however,  is  the  fact  that  various  optic  fow  vectors  could 
all  have  this  same  component.  The  second  component  of  the  optic  fow  (i.e.  the  one  that  is  perpendicular 
to  V7)  cannot  be  solved  for  using  only  the  intensity  constraint  equation.  At  least  one  more  equation  is 
needed  to  be  able  to  ..nd  the  complete  optic  X ow  v  as  a  unique  solution.  T  his  condition  is  what  makes  the 
problem  ill-posed  as  it  now  stands. 

There  are  four  principal  approaches  to  making  this  problem  one  that  can  be  considered  well-posed 
mathematically:  regularization,  correlation/ matching  methods,  spatiotemporal  ..Itering  methods,  and  ..[- 
tered  dinerential  methods  [3],  [2],  Sections  6  and  7  of  this  paper  examine  the  attempts  made  by  Bernard  and 
H orn  &  Schunck,  respectively.  Bernard's  approach  can  be  characterized  as  a  projected  dinerential  method 
which  employs  a  multiscale  technique.  According  to  this  method,  he  "hits"  the  governing  equation  with  a 
wavelet  basis,  makes  an  assumption  about  the  optic  fow  on  a  local  scale,  and  so  transforms  the  problem  into 
a  system  of  matrix  equations.  Thus  he  provides  a  satisfactory  number  of  equivalence  relations  to  uniquely 
solve  for  both  components  of  the  velocity  . . eld  [3],  [4],  [5],  [6],  [7],  In  a  completely  dinerent  approach,  H or n 
and  Schunck  rely  on  the  addition  of  a  smoothness  constraint  and  the  technique  of  calculus  of  variations  to 
regularize  their  problem  and  so  make  it  well-posed  [14], 


5.4  Additional  Issues  to  Consider 


Despite  the  existence  of  several  valid  methods  for  making  optic  fow  a  stable  and  solvable  problem, 
there  remain  a  few  additional  sources  of  error  which  require  addressing  [27,  pg.  70],  FI  ere  they  will  only  be 
mentioned  in  brief. 

The  ..rst  source  can  be  thought  of  in  terms  of  stochastic  error.  SpecL.cally,  the  digital  mapping  of  real- 
life  scenes  to  grayscale  image  sequences  will  invariably  include  the  presence  of  sensor  noise.  T  here  exists  a 
disparity  between  the  actual  velocity  of  the  objects  being  ..Imed  and  what  the  corresponding  velocity  will 
be  in  the  optic  fow  ..eld.  T  his  dinerence  may  be  described  as  a  random  noise  variable  which  is  statistical 
in  nature.  Various  techniques  exist  to  minimize  the  error  stemming  from  this  source  [23]. 

Another  key  issue  which  surfaces  is  the  fact  that  the  governing  equation  for  optic  fow  neglects  any 
interaction  by  second  order  terms  in  the  displacement  of  grayscale  intensity.  I  n  other  words,  it  only  accounts 
for  X  ow  velocity  and  pays  no  attention  to  the  enects  of  acceleration  within  the  X  ow.  T  his  leads  to  a  systematic 
error  whenever  the  magnitude  of  the  velocity  is  large  over  a  localized  area  [27,  pg.  70], 
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Finally,  we  note  that  error  could  arise  from  a  X  aw  in  the  model  for  optic  fow  itself.  Perhaps  the  physics 
and  the  theory  behind  the  governing  equation  don't  hold  and  the  assumption  that  the  brightness  of  an  object 
doesn't  change  over  time  (i.e.  that  grayscale  intensity  is  conserved)  is  not  valid.  These  problems  could  very 
well  manifest  themselves  in  a  negative  sense  as  one  goes  about  solving  for  optic  fow  depending  upon  the 
particular  image  sequence  [27,  pg.  70], 

Within  the  scope  of  this  paper,  methods  for  overcoming  the  ..rst  and  last  of  these  sources  of  error 
will  not  be  addressed  in  detail.  However,  the  second  issue  concerning  the  systematic  error  due  to  large 
displacements  is  dealt  with  under  B  ernard's  approach  in  the  next  section. 


6  Wavelet  Based  Multiscale  Approach 

6.1  W  hat  M  akes  Wavelets  U  nique 

As  was  stated  in  the  previous  section,  Bernard  attempts  to  overcome  the  ill-posedness  of  the  problem 
by  creating  a  system  of  matrix  equations  through  the  use  of  a  wavelet  basis.  But  what  makes  wavelets 
so  special  that  they  provide  a  better  way  to  solve  for  the  optic  fow  vectors?  What  is  special  about  the 
information  wavelets  can  extract? 

Bernard  describes  the  utility  of  wavelets  in  the  context  of  optic  fow  in  terms  of  their  ability  to  break 
through  the  challenge  posed  by  the  aperture  vs.  timealiasing  problem.  W  hen  trying  to  capture  the  motion 
of  an  object  in  an  image  sequence,  theditf  cuity  lies  in  ..nding  a  mathematical  "window"  which  is  big  enough 
to  measure  the  displacement  of  that  object  between  frames,  and  yet  not  so  large  that  the  motion  becomes 
so  i nsigni . .cant  that  it  is  undetectable  [3], 

An  illustrative  example  of  this  concept  presents  itself  if  one  thinks  of  watching  a  game  of  baseball 
under  two  dinerent  circumstances.  The  ..rst  involves  looking  at  the  ..eld  through  a  pair  of  binoculars. 
Although  one  is  able  to  zoom  in  and  watch  what  sign  the  catcher  is  giving  to  the  pitcher,  it  is  impossible 
to  simultaneously  see  that  the  runner  on  ..rst  base  is  moving  to  steal  second.  This  inability  to  see  large 
scale  motion  while  looking  through  a  small,  high  resolution  window  is  referred  to  as  the  aperture  problem. 
In  wect,  the  aperture  one  is  looking  through  is  not  of  su <t  dent  size  to  see  the  big  picture.  T he  second 
case  entails  watching  the  game  using  only  one's  normal  vision.  Now  the  viewer  can  see  that  the  runner  is 
attempting  to  steal  the  base,  however,  they  cannot  at  the  same  time  see  what  motion  the  catcher  makes  with 
his  ..ngers  as  a  sign  to  the  pitcher.  Here  the  window  is  so  large  that  small  scale  motion  becomes  i  nsigni.. cant 
and  so  is  lost  to  the  sensors  of  visual  perception.  T  his  is  known  as  time  aliasing  since  small  displacements 
become  aliased  out  of  the  picture  by  the  comparatively  large  and  low  resolution  viewing  window. 

In  terms  of  solving  for  optic  fow,  the  greatest  impediment  comes  when  large  scale  motion  is  occurring 
in  one  region  of  the  image  sequence  while  small  scale  motion  is  simultaneously  occurring  in  another.  T  he 
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following  discussion  will  demonstrate  the  properties  of  wavelets  that  allow  them  to  overcome  this  problem 
and  so  avoid  the  limitations  of  time  aliasing  vs.  aperture.  Upon  completion  of  this  general  overview  of 
wavelets,  we  will  see  how  Bernard  applies  them  to  the  speci..c  challenges  posed  by  optic  tow. 

6.2  Wavelets  and  M  ultiresolution  A  nalysis 

As  the  inherent  value  of  wavelets  stems  from  their  ability  to  perform  a  multiresolution  analysis  (M  RA), 
we  shall  start  there  by  de..ning  exactly  what  this  concept  means.  T  he  basic  principle  behind  a  multi  resolution 
analysis  is  to  decompose  a  function  space  into  a  sequence  of  subspaces,  denoted  by  Vj.  It  must  be  noted  that 
the  set  of  functions  mentioned  here  are  all  those  which  are  square  integrable,  namely  all  functions  /  such 
that  J  (f(t))2dt<  oo.  T  his  function  space  is  denoted  by  L2(K).  W  hen  referring  to  the  decomposition 
of  /,  it  may  be  helpful  to  think  of  the  function  as  a  signal  which  is  being  broken  down  into  its  components 
occurring  at  dinerent  scales.  The  subspace  V3  is  de..ned  according  to  its  integer  valued  resolution  j  such 
that  all  details  of  the  signal  which  are  of  scale  less  than  2'  j  are  suppressed  [19],  [28],  [26],  [16],  [20], 

A  number  of  requirements  are  placed  upon  the  subspaces  Vj  in  order  for  them  to  form  a  true  multireso¬ 
lution  analysis.  T  he  ..rst  states  that  V}  be  contained  in  every  subspace  which  is  higher  in  resolution  than  j. 
If  we  de..ne  the  jth  level  approximation  to  f(t)  as  fj(t )  then  /, (t)  e  Vj,  which  in  plain  language  means  that 
fj(t)  belongs  to  the  subspace  Vj.  Logically,  wesee  that  any  detailswhich  can  becaptured  at  resolution  level 
j  must  also  be  included  with  the  information  at  a  higher  resolution.  Therefore,  V:j  has  to  be  contained  in 
Vj+ 1,  which  is  represented  mathematically  by  V:j  c  Vj+i  for  all  integer  values  of  j.  Extending  this  principle, 
we  arrive  at  the  following  nesting  of  subspaces:  ...  c  Vfl  2  c  Vjl  i  c  Vj  c  V3+1  cVj+2c  ...  [22],  [28],  [19], 
[16]. 

The  details  which  are  lost  in  going  from  the  approximation  of  /  at  level  j  +  1  to  level  j  are  those  which 
occur  at  scale  2'  U'+D,  Wede..nethem  by  dj(t)  -  fj+i(t)  -  fj(t).  From  this  relationship  we  arrive  at  the 
equivalence  relation 

=  fj(t)  +dj(t). 

T  he  corresponding  subspaces  can  be  represented  accordingly, 

Vj+i  «  Vj  ®  Wj. 

Here  the  notation  Wj  refers  to  thedetail  space  at  resolution  level  j  and  Wj  is  orthogonal  to  Vj.  T  his  means 
that  if  we  took  the  inner  product  between  any  element  g  in  Wj  and  any  element  h  in  Vj,  the  result  would 
equal  0: 

h(t)g(t)dt~  0. 

T  his  breakdown  of  subspaces  can  continue  on  so  that  we  obtain  Equation  (7) , 

Vj+i  =  Wj  ®  Vj  =  Wj  ®  Wj ;  j  ©  Vj-,  i  =  ...  =  Wj  ©  Wjt  i  ©  Wj i  2  ©  -.  ©  Wj i  j  ©  Vj i  j.  (7) 
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If  Wj  is  orthogonal  to  V3  it  is  also  orthogonal  to  any  subspace  of  V3,  including  all  detail  spaces  Wk  such  that 
k  <  j.  From  Equation  (7),  we  see  that  our  approximation  space  at  level  j,  V3,  can  be  written  as  a  sum  of 
subspaces.  T  his  is  the  ..rst  requirement  for  an  M  RA  [28],  [19], 

The  second  necessary  condition  is  that  all  square  integrable  functions  can  be  represented  at  the  ..nest 
level  approximation  space  and  also  that  only  the  zero  function  is  contained  at  the  coarsest  level  of  resolution. 
Explanation  for  the  second  of  these  points  stems  from  the  fact  that  as  the  resolution  becomes  coarser  and 
coarser,  more  and  more  details  are  lost  and  at  the  limit  j  -s>=  -oc,  only  a  constant  function  remains.  In  order 
to  meet  the  requirement  of  square  integrability,  that  function  must  beexactly  equal  to  zero.  M  athematically 

this  can  be  represented  either  by  the  limit:  lim  V  =  {0}  or  by  the  in.mite  intersection  of  subspaces: 

o'-  i  i 

fi jVj  =  {0}.  Conversely,  as  resolution  increases,  more  and  more  details  are  included.  Therefore,  as 

the  resolution  level  approaches  positive  in . .n ity,  the  entire  space  of  square  integrable  functions  should  be 

recovered.  This  can  be  written  mathematically  either  by  the  limit:  lim  V3  =  L2(K)  or  by  the  closure  of 

_  o'-  1 

the  union  of  all  subspaces:  U, V7  =  L2(R)  [28],  [22], 

Thethird  requirement  states  that  all  the  spaces  {V3}  arescaled  versions  of  the  space  V0  which  is  named 

the  central  space.  If  we  say  that  the  function  f(t)  e  V3,  in  other  words  there  are  no  details  in  /(f)  which 

appear  at  scales  smaller  than  2'  j,  then  /(2t)  is  the  function  attained  by  compressing  f(t)  by  a  factor  of  2 

so  that  it  contains  no  ddtails  at  scales  less  than  2'  (j+1).  T  herefore,  /(2t)  e  V3+i-  This  condition  is  known 

as  scale  or  dilation  invariance  [28], 

Condition  number  four  mandates  that  the  space  V3  be  translation  or  shift  invariant.  Therefore,  if 
/(t)  e  Vo  then  also  f(t-  k)  e  Vo  for  all  integers  k.  The  combination  of  this  requirement  with  the  third  one 
relating  to  scale  invariance  give  way  to  the  fact  that  if  f(t)  e  Vo,  then  f{2jt-  k)  e  V3  [19],  [28],  [16], 

The  ..nal  requirement  for  a  multiresolution  analysis  is  that  there  must  exist  a  function  <j>  such  that  all  of 
its  translates  form  an  orthonormal  basis  for  Vo.  I  n  this  case,  orthonormality  requires  two  separate  conditions 
be^met,  one  of  orthogonality  and  theother  of  normalization  to  1.  T  he  ..rst  is  that  for  any  two  integers  k  ±  l, 
J  <j>(t  -  k)4>{t  -  l)dt  =  0  while  the  second  states  that  for  any  integer  k,  J  (</(f  -  k))2dt  =  1.  If  we 
apply  the  scale  invariance  property,  we  ..nd  that  {/(2t  -  k)}k  becomes  an  orthogonal  basis  for  Vl.  Using 
the  same  line  of  reasoning,  if  we  de..ne  <j)jk(t)  =  2-'/2</(2 jt  -  k),  then  {<frjk(t)}k  forms  an  orthonormal  basis 
for  the  space  V3.  In  referring  to  <p3k  as  a  basis  function,  we  mean  that  any  function  f3  which  lives  in  the 
space  Vj  can  be  represented  as  a  sum  of  these  basis  functions, 

fj(t)=  E  ajk<pjk(t). 

fc=i  l 

Extending  this  result,  we  say  that  any  square  integrable  function  /  which  necessarily  lives  in  L2(R)  can  be 
expanded  as  a  sum  of  the  basis  functions  across  all  scales, 

/W  =  E  E  ajk<t>jk  (t)- 

k=\  i  j=  i  i 

The  name  given  to  the  function  <j>  which  generates  the  basis  functions  for  all  spaces  {Vj}  is  the  scaling 
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function  of  the  multi  resolution  analysis  [19],  [28],  [16],  [26], 

In  order  to  construct  the  basis  functions  for  all  spaces  {Vj}  from  the  scaling  function  4>,  we  look  to 
an  equivalence  relation  known  as  the  Scali ng  Equation.  Its  development  begins  with  the  following  line  of 
reasoning.  According  to  our  previous  discussion  of  the  nesting  of  subspaces,  V0  c  Vi  and  therefore  any 
function  which  resides  in  V0  can  be  expanded  in  terms  of  the  basis  functions  for  Vi.  M  ore  speci . .cal ly,  the 
scaling  function  <f>eV o  can  be  written  in  terms  of  {<Plk(t)}k, 

i  i 

<t>(t)  =  X]  hk<t>ik{t)  =  ^2  ^2  hk<j)(2t-k).  (8) 

k=  i  1  fc=i  1 

Equation  (8)  is  called  the  Scaling  Equation.  T  his  same  process  can  be  applied  to  relatethe  scaling  functions 
between  any  two  successive  levels  of  resolution, 

<t>j\  i,m(f)  =  vV'  1<A(2Ji  h-m)  -  vV  J2  hk<j)(2Jt  -  k)  =  Y  hk<t>Jtk{t). 

k=\  1  k=  i  1 

The  ..Iter  coep  cients  {hk}  can  be  found,  using  the  fact  that  {4>ik(t)}  are  orthonormal,  by  taking  the 
inner  product 

hk  =  ifc)  =  VzJ  (j){t)4>{2t  -  k)dt. 

Using  this  process,  we  are  able  to  ..nd  the  functions  <j>jk{t)  which  form  a  basis  for  any  approximation  space 
Vj  [28],  [19],  [16], 

We  now  shift  our  attention  to  the  detail  spaces  denoted  by  {Wj}.  Ifweusethefactthat  Mm  V3  =  { 0} 
and  we  extend  Equation  (7)  to  lim  ,  we  ..nd  that  we  can  represent  the  approximation  space  Vj+i  by  the 

j'-  i  i 

orthogonal  sum 

Vj+i  =  ®i=i  i  Wk. 

If  we  let  j  ->■  oo,  we  see  that 

L2{ R)  =  ®U|  i  wk. 

In  other  words,  the  set  of  all  square  integrable  functions  can  be  broken  down  into  orthogonal  subspaces,  each 
of  which  can  detect  details  of  the  functions  up  to  a  given  resolution.  This  means  that  the  union  of  all  the 
basis  functions  for  {Wj}  serve  in  the  same  manner  as  a  basis  for  L2(K).  J  ust  as  the  approximation  spaces 
V:i  each  had  a  basis  {<frjk(t)}k,  every  detail  space  W3  in  a  multiresolution  analysis  has  its  own  orthonormal 
basis  {ipjk{t)}k.  In  similar  fashion  to  <pjk(t),  wede.me  ipjk{t)  =  2>/2ip(2H  -  k).  Thus  we  may  say  L2(R) 
has  an  orthonormal  basis  {tpjk(t)}jk  known  as  the  wavelet  basis,  also  referred  to  as  a  wavelet  family  [28], 
[19],  [16],  [20], 

The  function  4>{t)  which  is  used  to  construct  all  the  basis  functions  ikjk{t)  for  the  W3  spaces  is  known 
as  the  mother  wavelet.  By  de..nition  of  tpjk(t)  we  see  that  {V>( t  -  k)}k  belongs  to  W0  and  since  Equation 
(7)  tells  us  that  W0  c  Vi,  ip{t)  can  be  written  as  the  sum  of  the  basis  functions  for  Vi,  {4>lk{t)}k  : 
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l  _  l 

Ht)  =  Y  9  k4> (t)  —  V ^2  ^  ^  gk4>(2t  A).  (9) 

fc=i  l  fc=,  l 

Equation  (9)  is  known  as  the  W  avelet  Equation  and  it  can  be  extended  to  ..nd  the  relationship  between  the 
wavelet  and  the  scaling  function  at  the  next  ..ner  scale  for  any  two  successive  levels  of  resolution. 

J  ust  as  we  did  to  solve  for  the  ..Iter  coep  cients  of  the  scaling  equation,  we  apply  thefact  that  {4>lk{t)}k 
are  orthonormal  to  ..nd  the  ..Iter  coep  cients  gk, 

gk  =  (i>,<j>ik)  =  v7^  J  ip(t)(f>(2t-  k)dt. 

Using  this  process,  we  are  able  to  ..nd  all  of  the  functions  ipjk(t)  needed  to  form  a  basis  for  L2{ R)  [19],  [28], 
[16]. 

The  following  example  will  demonstrate  how  a  sample  square  integrable  function  can  be  decomposed 
into  a  set  of  approximations  and  details  occurring  at  various  levels  of  resolution  according  to  an  M  RA .  Our 
choice  of  scaling  function  and  wavelet  bases  is  referred  to  as  Daubechies  4  in  [26].  As  the  speci..cs  on  how 
to  construct  these  bases  are  described  by  Walker,  who  founds  his  work  upon  the  seminal  paper  by  Ingrid 

Daubechies  [8],  we  leave  the  reader  to  examine  these  references  so  as  not  to  draw  the  discussion  too  far 

(r) 

from  the  focus  of  how  a  multiresolution  analysis  works.  The  MATLABW  code  written  to  generate  these 
bases  and  the  ensuing  M  RA  of  the  samplesignal  is  listed  in  the  Appendices  of  this  paper  under  the  ..lename 
DaubA_  10. m  .  For  now,  let  it  sup  ce  to  say  that  these  scaling  functions  and  wavelets  were  designed  by 
Daubechies  to  be  orthonormal  and  of  compact  support  (i.e.  the  functions  only  have  non-zero  value  for  a 
..nite  closed  interval). 

(r) 

Figure  7  shows  a  sample  signal  which  is  represented  in  M  AT  LA  B  as  a  vector  of  length  1024. 


Sample  Signal  f(t)  =  20t2(1  -  t)4cos(12n  t) 


Figure  7:  Sample  Signal 
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Since  the  signal  /(f)  is  of  length  1024  =  210,  the  ..nest  resolution  level  to  which  we  can  decompose  it  is 
j  =  10  and  the  coarsest  I  a/ el  is  j  =  1.  The  coarse  to  ..ne  approximations  to  /(f),  /i(t)  through  /i0(f),  are 
displayed  in  Figures  8  and  9. 


Approximations  to  f(t)  at  Resolution  Level  j 


Figure  8:  Coarser  A pproximations  to  Sample  Signal 


Approximations  to  f(t)  at  Resolution  Level  j 


F  igure  9:  F  iner  A  pproximations  to  Sample  Signal 
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F  igures  10  and  11  show  the  details  of  f(t),  di(t)  to  dw(t),  from  coarse  to  ..ne. 


Details  of  f(t)  at  Resolution  Level  j 
Using  Daubechies-4  Wavelets 


0  0.25  0.5  0.75  1 


F  igure  10:  Coarser  Details  of  Sample  Signal 


Details  of  f(t)  at  Resolution  Level  j 
Using  Daubechies-4  Wavelets 


Figure  11:  F  iner  Details  of  Sample  Signal 
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From  examination  of  Figures  8  and  9,  one  can  see  that  the  ..nest  approximation  of  /(f)  at  resolution 
level  j  =  10  does  a  perfect  job  of  capturing  every  detail  about  the  signal.  This  corresponds  to  the  ..nest 
level  details  of  the  sample  signal,  dw{t),  taking  on  the  value  of  zero.  In  enect,  no  information  is  lost  when 
/(f)  is  approximated  by  /io(t).  As  the  approximations  get  progressively  coarser,  more  and  more  details  are 
left  out  until  we  see  that  the  coarsest  approximation  at  level  j  =  1,  only  a  constant  function  survives. 

We  now  seek  to  validate  the  assertion  of  Equation  (7)  and  prove  that  a  signal  can  be  perfectly  recon¬ 
structed  by  summing  up  the  approximation  to  it  at  level  j  -  J  with  the  details  from  level  j  down  to  j  -  J. 
We  choose  the  arbitrary  value  of  J  =  5  to  show  a  reconstruction  after  six  levels  of  decomposition.  Figure 
12  demonstrates  the  fact  that 

/(f)  =  dio(f)  +  dg(t )  +  d8(f)  +  (77(f)  +  de(t)  +  d5(t)  +  /5(f). 


Reconstruction  of  Sample  signal  f(t)  =  d10(t)+d9(t)+d8(t)+d7(t)+d6(t)+d5(t)+f5(t) 


F  igure  12:  Reconstruction  of  Sample  Signal 

So  far,  our  discussion  of  wavelets  and  multiresolution  analysis  has  taken  place  in  the  context  of  one¬ 
dimensional  function  spaces.  The  problem  of  optic  fow,  however,  involves  two-dimensional  images.  There¬ 
fore,  we  must  see  how  these  concepts  evolve  in  2-D. 

To  construct  a  two-dimensional  scaling  function,  one  merely  takes  the  tensor  product  of  two  one¬ 
dimensional  scaling  functions, 
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<j>(  XI, X2)  =  4>{xi)(t>(x2). 

In  turn,  the  Scaling  Equation  becomes 


1  1 

4>{xi,  X2)  =  2  ^2  ^2  hkj(j)(2xi  —  k,  2x2  —  l)  ■  (10) 

fc=i  1  /=  i  1 

Both  0(cci)  and  cf>(x2 )  satisfy  the  1-D  Scaling  Equation  <p(x)  =  \[2  J2  hk(j>(2x  -  k )  so  therefore  we  may 

it=i  1 

write  hk,i  =  hkhi.  T  hus,  the  2-D  Scaling  Equation  is  the  product  of  two  1-D  Scaling  Equations  [28], 

We  now  approach  a  shift  in  terminology  for  the  two-dimensional  case.  Here,  a  2-D  M  other  Wavelet  can 
beany  of  four  combinations  denoted  by, 

01(x1,x2)  =  tp(xi)4>{x2), 

02(xi,x2)  =  <p(xi)ip(x2), 

03(x  !,x2)  =  i)(x\)ij)(x2), 

04(x  l,x2)  =  4>(xi)(j)(x2). 

Notice  that  the  general  2-D  Mother  Wavelet  0n(xi,x2 ),  with  index  n  =  1,2,3,  or  4,  can  be  composed  of  a 
tensor  product  of  either  1-D  scaling  functions,  wavelets,  or  both.  T  hesigni..cance  of  each  lies  in  the  dinerent 
types  of  information  they  provide.  61  is  known  as  the  High-Low  component  because  it  provides  detailed 
high-resolution  information  in  the  x1  direction  and  long-term  trend  information  in  the  x2  direction.  This 
stems  from  the  fact  that  1-D  wavelets  form  bases  for  the  detail  spaces  of  a  signal's  decomposition  while  1-D 
scaling  functions  form  bases  for  the  approximation  or  "trend"  spaces.  Thinking  of  an  image  as  a  2-D  signal, 
61  would  serve  best  to  capture  vertical  lines  appearing  in  the  image.  Similarly  62,  being  the  Low-High 
component,  would  provide  the  most  information  about  horizontal  lines  occurring  in  an  image  [3,  C  h.  4],  [28], 

At  this  point,  the  basic  theory  behind  waveldts  and  muitiresolution  analysis  has  been  explained.  We 
now  seek  to  demonstrate  how  B  ernard  makes  use  of  these  concepts  to  solve  the  problem  of  optic  fow. 

6.3  Application  of  Wavelets  to  Optic  Flow 

A  key  question  still  remains  unanswered.  What  is  the  exact  purpose  wavelets  serve  in  solving  for  the 
coep  dents  of  the  optic  fow  vectors  v\  and  v2l  The  answer  lies  in  the  multiscale  nature  of  wavelets. 

In  real  life,  actual  motion  of  an  object  happens  over  a  continuum  of  both  space  and  time.  This  is  not 
thecasein  thecontext  of  an  image  sequence.  In  thisdigital  environment,  motion  is  represented  by  an  object 
comprising  a  set  number  of  pixels  that  shifts  over  a  certain  number  of  pixels  in  the  xi  and  x2  directions.  It 
may  help  to  think  of  the  space-time  continuum  being  projected  onto  a  discrete  grid  of  pixels  and  frames. 
T  his  is  exactly  why  wavelets  become  so  bene..cial. 

We  de..ne  the  2-D  wavelet  family 

#"k(x)  =  2 i/26n  (2JX  —  k)  where  j  £  Z  and  k  e  Z2  and  n  =  1, ...,  4. 
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Here  Z  represents  the  set  of  all  integers  and  Z2  the  2-D  Cartesian  grid  of  integers.  Let  us  consider  a  video 
sequence  where  every  image  is  of  a  2J  by  2J  pixel  grid.  Standard  practice  is  to  think  of  this  area  being 
normalized  so  that  it  has  length  1  in  both  the  x\  and  the  x2  direction.  Accordingly,  each  pixel  is  considered 
to  be  a  square  with  each  of  its  sides  having  length  =  2'  J.  At  each  resolution  level  j,  the  2-D  wavelet 
cannot  detect  details  of  the  image  which  are  of  scale  less  than  2'  j.  In  other  words,  0"k(x)  can  only  detect 
features  of  an  image  which  are  of  size  2Ji  j  pixels  by  2Ji  j  pixels  or  larger.  The  wavelets  at  the  ..nest 
resolution  level  J,  namely  0"k(x),  can  therefore  zoom  in  on  details  which  are  only  one  pixel  wide  by  one 
pixel  long  [3],  [7], 

From  this  result  we  are  able  to  see  how  a  multiresolution  analysis  performed  with  a  two-dimensional 
wavelet  basis  allows  us  to  capture  both  large  and  small  scale  motion  occurring  in  an  image  sequence  at  the 
same  time.  An  object  which  is  moving  at  a  relatively  high  velocity  will  be  represented  as  jumping  across 
large  gaps  in  the  pixel  grid  from  frame  to  frame.  Wavelets  at  the  coarser  levels  of  resolution  have  a  larger 
support  and  so  they  will  be  able  to  detect  this  change  in  location  of  the  object.  Conversely,  wavelets  at 
..ner  levels  of  resolution  do  not  have  a  support  large  enough  to  enclose  the  entire  displacement  travelled  by 
the  object  and  so  they  cannot  measure  its  motion  within  their  window.  This  is  the  aperture  problem.  0  n 
the  other  hand,  an  object  travelling  at  a  relatively  low  velocity  will  only  move  across  a  few  pixels  in  the 
grid  from  frame  to  frame.  Wavelets  at  the  ..ner  levels  of  resolution  will  have  support  small  enough  to  zoom 
in  and  detect  the  fact  that  the  object  may  have  only  shifted  over  by  a  few  pixels.  For  those  wavelets  at 
the  coarser  I  e/els  of  resolution,  however,  the  motion  may  not  represent  a  large  enough  trend  and  so  it  will 
suppressed  as  too  small  of  a  detail.  T  his  is  what  is  known  as  time  aliasing. 

We  see,  however,  that  the  time  aliasing  vs.  aperture  problem  is  avoided  when  we  simultaneously  apply 
all  levels  of  resolution  in  our  M  RA  to  the  motion  occurring  between  two  frames.  The  large  scale  motion 
is  represented  in  the  coarser  levels  of  resolution.  This  information  will  translate  into  long  vector  arrows 
in  the  optic  fow  ..eld.  Small  scale  motion  can  be  seen  at  the  ..ner  levels  of  resolution  and  will  translate 
into  shorter  vector  arrows  in  the  solution  for  the  optic  fow.  Additionally,  the  application  of  the  dinerent 
components  of  the  2-D  wavelets,  namely  6>1,02,6>3,  and  94,  to  the  image  sequence  allows  us  to  account  for  the 
detection  of  an  object's  motion  in  any  direction  [3],  [4],  [5],  [6],  [7], 

In  theory,  everything  seems  to  work  out.  Let  us  now  turn  to  Bernard's  method  of  solving  for  the 
optic  fow  vectors  to  see  if  these  wavelets  really  are  overcoming  the  challenges  of  the  problem  when  put  into 
practice.  Keep  in  mind  throughout  the  discussion  that  our  end  goal  is  to  solve  for  a  ..nite  number  of  vectors, 
each  of  which  will  be  applied  to  a  certain  subgrid  of  pixels  in  the  image  to  represent  the  optic  fow. 

Once  again  we  begin  with  the  governing  equation  for  optic  fow, 


dl  .  .  dl  ,  dl 

—„l(x)+— «(><)+ -  =  °. 


We  hit  the  governing  equation  with  a  basis  function  =  6n(x  -  u)  de..ned  in  our  notation  earlier,  and 
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then  integrate, 

II  Idx^Vl^  +  'dIV2^  +  %S)0n^X~U^dXldX2  =  0’  for  n  =  l,..,N  .  (11) 

Here  TV  could  equal  3  or  4.  Equation  (11)  can  be  expanded  and  rewritten  in  inner  product  notation, 

(dk"h^)  +  (akV2X)  +  ( I'"1)'0’  torn=l,..,N  .  (12| 

At  this  point,  Bernard  makes  his  most  critical  assumption.  He  states  that  the  optic  tow  vectors  vi(x) 
and  V2 (x)  are  constant  over  the  supports  of  the  basis  functions  0™  ,  i.e. 

i!i(x)  =  wi(u)  for  all  X  e  support  9™,  for  all  n , 

V2(x )  =  i>2(u)  for  all  X  e  support  9™,  for  all  n. 

Therefore,  if  the  optic  tow  vectors  are  constant  over  the  domain  of  the  inner  products,  they  can  be  taken 

outside  of  the  integrals.  In  addition,  since  the  inner  products  are  integrals  being  taken  with  respect  to  x, 

d 

we  can  take  the  partial  derivative  with  respect  to  time  —  outside  of  the  third  inner  product  term, 

at 

+  |  M)  -  °,  f  or  n  =  1 , . N  .  (13) 

If  we  perform  an  integration  by  parts  for  the  ..rst  two  inner  products,  we  arrive  at  the  equation 


for  n  =  1, TV 


Note:  The  terms  16™  are  equal  to  0  and  therefore,  they  do  not  appear  in  Equation  (14).  This 

b 

is  because  a  and  b  represent  the  boundaries  of  the  image  and  the  basis  functions  6™  only  have  ..nite  support 
so  they  vanish  on  the  boundaries. 

Equation  (14)  is  a  key  equation  because  if  we  now  turn  back  to  9 ^  for  n  =  1  ,..,N  where  TV  =  3  or  4  as 
de..ned  previously  in  our  notation,  we  have  a  projected  system  of  TV  equations  with  two  unknowns  vi(u) 
and  v2(u), 
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We  say  that  this  system,  which  Bernard  refers  to  as  ( S ),  is  projected  because  for  each  u  that  we  could 
choose  in  translating  across  the  image,  we  are  able  to  say  that  the  optic  fow  vector  v  is  constant  over  the 
support  of  6 "  [3],  The  problem  with  this,  however,  is  that  u  lives  in  the  continuum  so  therefore  our 

individual  optic  fow  vectors  could  be  constant  over  a  certain  area  of  the  image  but  not  necessarily  for  a 
certain  group  of  pixels.  T  his  leaves  us  with  the  fact  that  westill  have  a  vector  ..eld  v  we  aresolving  for  and 
not  a  ..nite  system  of  individual  vectors  as  we  desire. 

Furthermore,  there  is  no  mention  of  scale  in  (S).  We  need  to  ..t  in  some  way  of  re.ming  this  scheme 
so  that  we  can  account  for  average  trends  and  dinerences  over  various  groupings  of  pixels.  This  insertion 
of  multiresolution  analysis  will  allow  our  discrete  set  of  optic  fow  vectors  more  closely  approximate  the 
continuous  optic  fow  ..eld  [3], 

Finally,  we  must  also  contend  with  the  issue  that  an  image  sequence  is  sampled  in  time  ( i.e.  motion  is 
captured  with  a  discrete  number  of  image  frames).  T  herefore,  the  right  hand  side  of  ( S )  which  is  a  partial 
derivative  with  respect  to  time, 

must  be  approximated  with  some  sort  of  ..nite  dinerence  such  as: 

dl 

—  ~i(t  +  i)-l(t), 

where  t  +  1  really  just  represents  the  next  image  frame  from  the  one  at  time  t. 

In  his  thesis  [3],  Bernard  actually  uses  a  higher  order  estimate  to  the  time  derivative  and  measures  the 
optic  fow  at  each  t  +  1/2  according  to 


dl{t+  1/2) 
dt 


(15) 


This  approximation  makes  sense  if  one  considers  the  Fundamental  Theorem  of  Calculus  for  derivatives: 
f°(c)(b  -  a)  =  f(b)  -  f(a)  where  b  —  a  =  1  for  our  example. 

Now  we  will  introduce  dinerent  scalings  of  the  basis  functions  so  that  we  may  perform  an  M  RA  later 
on  to  solve  for  our  optic  fow  vector: 

C(x)  =  *'  lr  (~~)  ’ 

where  s  is  a  continuous  scaling  index.  We  examine  the  inner  products  on  the  left  hand  side  of  (5)  evaluated 
at  t  +  1/2, 


(16) 
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We  use  the  following  averaging  estimation  to  approximate  I  at  a  non-integer  time  (i.e.  in  between  frames): 


lit  + 1/2) 


/( t)  j- 1 [t  4- 1) 
2 


(17) 


If  we  make  use  of  (15),  (16),  and  (17)  we  can  form  a  new  projected  system  of  equations  that  includes 
various  scales  of  the  basis  functions  and  also  is  discretized  in  time, 


\  '  /  7 ( /  +  1 )  +  l{t)  dflL 
,fr2  \  2  ’  dxi 

/Ut  +  p  +  Ut)  dC 
A  A  2  ’  dxt 


vi(u)=\l(t+l)-I(t),61u ) 


(775) 


One  problem  still  remains,  we  are  still  trying  to  solving  for  our  optic  fow  as  a  vector  ..eld  rather  than 
a  set  of  individual  vectors  [3],  [5],  [6],  [7],  This  is  because  in  (775)  ,  we  are  still  employing  the  use  of  a 
continuous  translation  parameter  u  and  a  continuous  scaling/  dilation  parameter  s  .  In  order  to  resolve  this 
issue,  we  sample  the  continuous  parameters  u  and  s  to  create  discrete  parameters  k  and  j  .  We  formulate  a 
new  set  of  basis  functions  which  Bernard  now  begins  calling  a  wavelet  family  {0"k}„=i,..,;v;  ji  z:k2z2  de..ned 
by: 


6»Jk(x)  =  Z>/20n(Z>x  -  k), 


where  /  is  a  resolution  index,  and  k  =  (fci,  fe2)  is  a  2-D  translational  index.  The2-D  wavelet  6>”k  is  located 
around  (2>  jfci,  2'  ffc2)  and  has  support  over  a  domain  of  size  proportional  to  2'  i  . 

Using  this  new  set  of  discrete  translation  indices,  we  may  now  express  our  system  of  equations  as: 


Ml 

dxi 


de% 

dx\ 


vi  +  (  I 


dO} k’ 

dX2 


vi  +  (  7, 


del 

dX2 


V2  = 


V2  = 


dA* 


—  eN 


(5j-k ) 


T  he  approximations  for  the  discretization  of  time  used  in  (775)  are  left  out  in  (5/k)  merely  to  reduce  clutter 
in  the  equations  and  allow  the  reader  to  see  how  we  now  have  a  ..nite  system  of  equations  (either  3  or  4) 
for  which  we  must  solve  for  2  variables  v\  and  w2  for  each  translational  and  resolutional  combination  j, k  . 
The  ..nite  dinerence  approximations  to  the  time  derivatives  used  in  (775)  are  actually  used  for  performing 
the  calculations. 
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As  such,  solving  for  the  optic  X ow  has  gone  from  being  an  ill-posed  problem  to  one  that  is  well-posed. 
T  hus,  we  now  have  a  matrix  manipulation  problem  to  solve  for  a  ..nite  number  of  vectors  since  j  and  k  are 
..nite  [3],  [7]: 


(18) 


Seemingly,  all  of  our  obstacles  have  been  overcome  in  transitioning  to  the  ..nite  system  of  dinerence 
equations  in  Equation  (18).  The  trick  to  seeing  that  the  solutions  to  the  problem  in  Equation  (18)  will 
be  stable  (i.e.  the  time-aliasing  vs.  aperture  problem  has  been  surmounted)  comes  from  the  fact  that  the 
velocity  vectors  v(x,t)  are  being  solved  for  at  each  resolution  level  j  in  the  multiresolution  analysis.  The 
determinant  of  could  be  0  or  very  close  to  it  if  the  motion  is  occurring  on  a  scale  much  dinerent  than 
2'  j  in  that  region  covered  by  the  support  of  9]k.  This  would  cause  the  approximation  to  v  to  be  very 
unstable  at  that  level  of  resolution.  However,  due  to  the  fact  that  Bernard  uses  two  indices,  j  and  k,  to 
achieve  a  sampling  which  is  spatially  localized  at  various  scales,  the  rate  of  approximation  for  a  given  region 
is  unhindered  by  a  particular  resolution  level  yielding  an  unstable  result.  In  other  words,  since  there  can 
only  be  one  vector  assigned  to  represent  the  optic  X ow  at  location  x  at  time  t,  the  solutions  from  each  level 
of  resolution  collectively  represent  enough  information  to  describe  whatever  motion  is  occurring,  either  large 
or  small  scale,  with  one  single  vector.  T  hus  we  see  that  Bernard's  approach  can  be  viewed  mathematically 
as  an  attempt  to  stabilize  the  approximation  to  v  more  so  than  as  a  way  to  regularize  the  problem  [3],  [7], 
Issues  of  stability  with  regards  to  solving  an  overdetermined  system,  as  Equation  (18)  plainly  is,  are  not 
addressed  within  the  scope  of  this  paper,  however,  they  are  dealt  with  briefy  by  Bernard  in  [3,  Ch.  4.1.4], 
The  ..nal  challenge  to  implementing  Bernard's  method  of  solving  for  optic  fow  stems  from  having  to  take 
the  derivative  of  a  two-dimensional  wavelet  as  is  required  to  set  up  the  system  of  equations  Sj k. 


6.4  Obstacle  to  I  mplementation 


Although  the  computation  of  the  inner  products  (  I, 


and  ( I 


38% 


may  not  seem  overly  compli- 


\"  ’  dx\  /  \  ’  qx 2 

cated  at  ..rst  glance,  a  signi ..cant  obstacle  lies  in  the  requirement  to  take  the  derivative  of  a  two-dimensional 
wavelet.  Here  the  di<t  culty  arises  from  the  need  to  construct  wavelets  in  a  manner  that  lends  itself  to 
taking  their  derivative.  Our  contention  is  that  we  need  not  construct  the  derivative  at  every  point  over  R2. 


Rather,  when  dealing  with  a  discrete  signal  I,  the  value  of 
points  of  (xi, X2)  [22], 


<9(9 


j  k 


dxi 


and 
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ik 


dxi 


need  only  be  known  at  the  integer 


The  discussion  of  this  theory  [22]  will  be  conducted  in  one-dimension  so  as  to  avoid  the  confusion  and 
clutter  accompanying  the  analogous  explanation  in  two-dimensions.  Our  example  problem  begins  with  g(t), 
a  1-D  analog  signal  which  has  compact  support  [0,  L\  where  L  =  2J  for  some  number  J  and  g(L)  =  0.  T  he 
discrete  approximation  to  g(t)  is  formed  so  that  the  sampling  interval  is  of  length  1,  b(n)  =  g(n).  To  ease  the 
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burden  of  notation,  we  employ  the  normalized  signal  /(f)  having  support  0  <  t  <  1  such  that  /(f)  =  g(2Jt). 

1  /  Tl  \ 

T  herefore,  we  see  that  /(f)  is  sampled  at  intervals  of  length  2>  J  =  —  so  that  /  ( —  J  =  b(n). 

Analysis  of  thesignal  will  be  conducted  using  a  dyadic,  multi-resolutional  decomposition  of  L2{ R)  based 
upon  a  scaling  function  <f>  of  compact  support  [0,1V]  with  /(IV)  =  0.  By  dyadic,  we  mean  that  each  level 
of  resolution  drops  down  in  coarseness  by  a  factor  of  2.  As  is  stated  in  [19],  /  satis.. es  the  scaling  equation 
built  from  a  discrete  ..iter  {hn}  having  support  0  <  n  <  N  -  1, 

_JVj  1 

/(t)  =  ^2  £  M(2t-  0- 
1=0 

In  the  multi-resolutional  decomposition  of  L2(R)  based  upon  4>,  thefunction  /(f)  is  said  to  belong  to  the 
space  Vj  where  2'  J  =  y  is  the  sampling  interval.  The  goal  of  our  approach  is  to  initialize  the  multi  resolution 

Jj 

{/  Q\  'j 

(/,  Uj  k)  )  \  which  are  indexed  over  a  range  min  <  k  <  max  of  non- 

\  'It  fc=min 

zero  values,  and  then  to  engage  a  cascade  algorithm  which  is  based  upon  the  scaling  equation  to  iteratively 
..nd  the  inner  products  j^/,  j)  kJ  for  1  <  j  <  J. 

From  [8],  the  scaling  equation  can  be  evaluated  for  an  integer  value  according  to 

_JVi  1 

(j>(n)  =  V2  hi4>(2n  -  l) 

1=0 

=  s/2  E  hki  2  n<t>(k)- 
k 

where  k  is  indexed  by  0  <  k  -=2n  <  N  -  1. 

De..ning  the  matrix 

mOn,k  =  | 

and  restricting  to  the  integer  values,  the  scaling  equation  can  be  rewritten 


s/ 2  hk  j  2n 
0 


0<fc  —  2n  <  N  —  1 
otherwise 


(j>{n)  =  E  ?n0 n,k<l>(k)  where  0  <  k  —  2n  <  N  —  1. 


(19) 


Dinerentiating  the  scaling  equation,  a  new  scaling  equation  for  results, 


N  i  1 


/°(t)  =  2V2  hi<j>\2t-  l). 


(20) 


1=0 


If  we  restrict  Equation  (20)  to  the  integers  and  apply  the  principle  of  substitution  from  Equation  (19)  we 
obtain 


E  </°(n)  =  ^2  ??i0n,fc/°(fc) 

'  b 


where  0  <  k  —  2n  <  N  —  1. 


(21) 


31 

From  Equations  ( 19)  and  (21),  we  see  that  the  sequence  of  integer-sampled  values  of  the  scaling  function 
< i>  comprise  a  right  eigenvector  for  the  matrix  mO  with  eigenvalue  1,  and  the  integer- sampled  values  of  the 
derivative  of  the  scaling  function  </>°  make  up  a  right  eigenvector  for  ??i0  with  eigenvalue  1/2.  In  [8],  a 

rather  cryptic  process  unfolds  in  order  to  ..nd  which  eigenvectors  correspond  to  the  sampling  of  the  scaling 

(r) 

function  and  the  sampling  of  its  derivative.  The  exact  code  required  to  perform  this  process  in  M  AT  LA 
is  included  in  the  Appendix.  In  short,  the  ..rst  part  of  Daubechies'  scheme  chooses  that  A  =  1  eigenvector 
which  is  orthogonal  to  the  left  eigenvector  [1, 1,1, 1]  of  m 0  whose  sum  of  components  is  equal  to  1.  This 
will  be  the  sequence  (</>(n)}  where  0  <  n  <  N  -  1  which  are  the  non-zero  values  of  the  integer  samplings 
of  the  scaling  function.  The  second  part  of  the  scheme  chooses  that  A  =  1/2  eigenvector  whose  sum  with 
the  vector  [1,2, 3, 4]  belongs  to  the  span  of  the  previous  left  eigenvector  [1,1, 1,1],  This  will  in  turn  be 
the  sequence  (</°(?r)}  where  0  <  n  <  N  -  1  which  are  the  non-zero  values  of  the  integer  samplings  of  the 
derivative  of  the  scaling  function  [8], 

Having  found  thesequence  {/><V)}„=o,..,jVi  i,  we  may  initialize  the  algorithm  for  computing  inner  prod¬ 
ucts  of  the  signal  with  the  derivative  of  the  scaling  function.  We  continue  to  use  J  as  the  dyadic  exponent 
determined  by  the  support  of  the  original  signal  g(t),  so  that  2 J  =  L.  In  addition,  we  make  use  of  the 
de.mition  ^°J>TO(f)  =  2J/2</°(2Jt  -  to).  For  any  integer  value  of  to,  we  can  approximate  the  value  of  the 
inner  product  ^/,  (</JiTn)°)  according  to: 

i 

(/>  (</j,  m) °)  =  J  fit )  --Jt  to))  eft, 

t=i  1 

using  the  Chain  Rule  and  the  fact  that  0  <  t  <  1, 

i 

=  V2J2J  J  f(t)(/)0(2Jt  —  m)clt , 

t=  0 


breaking  up  the  integral  into  segments, 


_  2J  \  i 

=  V2J2J 

n=  0 


(n  +  l)/2J 

J  f(t)<p%2Jt  -  m)dt, 

t=n/2J 


using  a  R  iemann  approximation  to  the  segments  of  the  integral, 

2Ji  1 

=*  V272J  Y  f  (n/2J)  •  (j>°{ 2J  (n/2J)  -  to)  •  (1/2J)  , 

n  =  0 


cancelling  out  the  2J  factors  and  substituting  /  (n  /2J)  =  bin), 

_ 2Ji  1 

\Fif  Y^  ^(n)  ■  <fXn— 


n  =  0 


to), 


rsj 
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making  use  of  the  convolution  operator  *, 

=  \p2J  ■  (h*  (ji  ^  (m), 

where 

h  =  {b(n)}nnZL0l  1  =  2'  1 

&  =  {*(")}„  =  0 
^°(n)  =  (jP(—n). 

Thus  the  inner  product  of  our  signal  with  the  derivative  of  the  scaling  function  has  been  reduced  to 
a  seemingly  simple  convolution  between  the  integer-sampling  points  of  these  two  functions.  However,  one 
issue  still  remains  before  we  can  begin  computing  these  coep  cients,  that  being  the  boundary  conditions  of 
the  signal.  Obviously,  a  translated  and  dilated  derivative  scaling  function  will  have  dinerent  boundaries 
than  the  original  <y>°.  This  means  that  their  support  will  extend  past  the  boundaries  of  the  signal,  drastically 
anecting  the  corresponding  inner  products,  a  fact  that  cannot  be  ignored  [22], 

The  basis  for  addressing  this  problem  is  presented  in  [19,  Ch.  7.5],  He  contends  that  in  a  case  such 
as  this,  one  must  employ  a  circular  convolution.  The  idea  is  to  extend  the  "discrete  signal"  h  to  make  it 
a  periodic  signal  with  period  L,  and  the  "discrete  ..Iter"  _^°  to  a  periodic  ..iter  with  period  ma x(N,L).  In 
essence  the  convolution  becomes  a  "circular  convolution"  of  two  periodic  signals  that  has  itself  a  fundamental 
period  L.  From  this  information,  we  extract  the  "fundamental  period  signal",  a  sequence  indexed  by 
0  <  k  <  N  —  1.  This  sequence  represents  what  we  were  seeking,  the  inner  products  at  resolution  level  J 
[24],  [22], 

Within  the  WAVELAB^  add-on  package  to  MATLAB  ,  written  by  Dr.  David  Donoho  and  his 
colleagues  at  Stanford  University,  there  exists  an  m-..Le  aconv.m  which  performs  the  circular  convolution 
according  to: 


J  Level  Inner  Products  =  V2J  ■  aconv{DphiO,b). 

where  DphiO  represent  the  integer  samplings  of  thederivative  of  the  scaling  function,  {<f>°(?r)}„=o,..,wi  i<  ar|d 
b  is  the  vector  containing  the  integer  samplings  of  the  original  function  g(t)  [22], 

Having  found  the  inner  products  at  level  ,J,  denoted  from  now  on  as  J  Levi  nnerP  rods,  we  seek  now  to 
..nd  the  inner  products  at  resolution  level  J-  1,  denoted  by  J  _  lLevi  nnerP  rods.  We  make  use  of  the  scaling 
equation  as  it  applies  to  any  two  successive  levels  in  a  multiresolution  [19], 

N\  1 

^Ui  1),  m(^)  =  bf,  2n0ji/(t), 

1=0 


(22) 


33 

evoking  once  more  the  de.  .n  it  ion  4>jk(t)  =  2^24>(2 H  -m).  We  differentiate  Equation  (22)  to  yield 

o  N'  1 

(</>(  J\  1),  m)  (*)  =  I]  H  2 n  {4>J,  l )  0  (t)  ■  (  23) 

1  =  0 

Notethat  thederivativein  Equation  (23)  isbeingtaken  after  thescaling function  $  isdilated  and  translated, 
i.e.  (4>j'k)0  not  (0°)  .  Substituting  this  result  into  the  inner  product  we  see 


Here,  let  h  denote  the  sequence  of  non-zero  coep  cients  for  the  ..Iter  h(k)  =  hk  indexed  over  the  range 
0  <  k  <  N  -  1.  Further,  let  A  denote  the  ..iter  under  involution,  i.e.  h(k)  =  A(-fc)  which  is  a  sequence 
indexed  over  -(N  - 1)  <  k  <  0.  We  apply  our  ..nding  of  J  LevinnerP  rods  to  rewrite  Equation  (24)  in  terms 
of  a  convolution  as  outlined  in  [19], 

J  _  lLevi  nnerProds  =  ^/,  {^(J.  1(j  ^  Levi nnerP rods  *h)j  (2m).  (25) 

Notice  that  after  the  convolution  is  performed,  the  result  is  downsampled  by  a  factor  of  2  according  to  the 
argument  2m.  This  dyadic  decimation  of  the  inner  products  corresponds  to  thecoarsening  of  the  resolution 
by  a  factor  of  2  as  one  goes  from  level  J  to  level  J  -  1. 

The  issue  of  the  signal's  boundary  must  also  be  remembered  when  computing  these  lower  resolution 

(r) 

inner  products.  Fortunately,  the  WAVELABw  m-..le  DownDyadLo.m  accounts  for  this  fact  in  addition 
to  performing  the  required  downsampling.  Therefore,  we  can  e<t  ciently  engage  our  cascade  algorithm  for 
..nding  the  J-l  level  inner  products  by 

J  _  lLevi  nnerP  rods  =  DawnDyadLoQ  Levi  nnerP  rods,  h),  (26) 

where  h  represents  the  vector  containing  the  values ij(fc)  =  hk  indexed  over  the  range  0  <  k  <  N  -  1. 

The  cascade  algorithm  based  upon  the  scaling  equation  extends  further  down  the  levels  of  resolution  as 

well, 

J  _  2Levlnner  Prods  =  ^/,  2)  TO)  ^  =  (j  _  lLevi  nnerP  rods  *h j  (2m),  (27) 

(?) 

and  is  implemented  in  WAVELAB  by 


J  _  2Levl  nner  P  rods  =  DownDyadLo(]  _  lLevi  nnerP  rods,  h). 


(28) 
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According  to  the  dyadic  decimation  of  the  inner  products,  this  process  of  performing  a  circular  convolution 
and  then  downsampling  allows  the  muitiresolution  to  extend  down  J  levels,  i.e.  to  the  j  -  j  =  o  level  of 
resolution,  where  there  will  only  be  1  inner  product  being  computed  for  the  set  J_J  Levi  nnerP  rods  [22], 
Throughout  this  discussion,  we  have  been  using  the  scaling  function  p  to  generate  a  basis  for  our 
multiresolution  analysis.  What  if  we  had  desired  to  use  the  mother  wavelet  ip  instead?  Looking  to  the 
wavelet  equation,  we  see  that  the  answer  is  not  overly  complicated.  Recall, 

_JVi  1 

ip(t)  mV?  J2  gip>(2t  —  l), 

1  =  0 

where  the  ..Iter  coep  dents  gt  can  be  solved  for  gt  =  -Jl  f  ip(t)<p(2t  -  l)dt. 

J  i  i 

We  restrict  the  wavelet  equation  to  the  integer  values  as  before  in  the  scaling  equation, 

,-N  i  1 

ip(n)  =  V2  J2  gi4>{2n-  l), 

1  =  0 


=  ^1®!  2 n<P{k), 

k 

where  0  <  k  -  2n  <  N  -  1.  Now  if  we  de..ne  the  matrix  ml 

_f  y/2gkl  2n  0  <  k  -  2n  <  N  -  1 
171  n’k  \  0  otherwise 

we  may  rewrite  the  wavelet  equation  restricted  to  the  integers  as 

4’{n)  =  2  ml n,k<t>(k)  where  0  <  k  —  2n  <  N  —  1. 
k 

Now,  by  applying  the  results  of  Equations  (19)  to  (21),  we  see  that  we  may  ..nd  the  exact  values  of  the 
derivative  of  the  mother  wavelet  at  the  integers  according  to  the  matrix  multiplication 

ip°(n)  =  2  J2  ml n,k<P°(k)  where  0  <  k  —  2n  <  N  —  1. 

k 

We  now  seek  to  ..nd  the  inner  products  of  the  signal  with  these  values  {^°(n)}„= i.  This  can  be 
accomplished  using  the  same  process  as  outlined  for  ..nding  the  J  level  inner  products  of  <p°.  T  herefore,  the 
circular  convolution  is  performed  to  ..nd 

Wavelet  J  Level  Inner  Products  =  ^/,  (V’j,  m)°)  =  ^2J  ■  aconv(DpsiO,b), 

where  DpsiO  represents  the  integer  samplings  of  the  derivative  of  the  mother  wavelet,  {ip°(n)}n=o:..,Ni  i,  and 
b  is  the  vector  containing  the  integer  samplings  of  the  original  function  g(t)  [22], 
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J  ust  as  in  Equation  (22),  we  can  extend  the  wavelet  equation  to  any  two  successive  levels  in  the  multireso¬ 
lution  [19], 


N  j  1 

V’Ui  i) ,m(*)  =  9l'<  (29) 

1  =  0 

again  using  the  de. .nition  ipjtk(t)  =  2j/2tp(2 jt- m)  [19].  We  apply  the  results  of  Equations  (23)  and  (24)  to 
obtain  the  equation  for  inner  products  with  wavelets  at  lower  levels  of  resolution, 

(/,  (^(Ji  1) ,  m)  )  =  M  9l\  2 n  (/,  (^J,z)°)  • 

'  '  1=0 

Since  ip  is  built  from  cp  as  in  Equation  (29),  the  exact  same  initialization  to  the  cascade  algorithm 
applies.  The  integer  samplings  of  the  derivative  of  the  scaling  function,  1,  are  used  to  ..nd  the 

inner  products  at  level  J,  //,  (</>7>  m)0)'  what  we  have  named  J  Levi nner Prods.  These  areused  by  the  m-..le 

DownDyadLo.m  to  compute  the  J  -  1  level  inner  products  with  the  scaling  function,  ^/,  ^(Ji  1(i  m) 

However,  WAVELAB®'s  matching  m-..le  DownDyadHi.m  uses  them  also  to  compute  the  J  -  1  level  inner 
products  with  wavelets,  (f,(ip(J.  1}  TO)  \.  Within  the  process  of  evaluating 


W  avj  _  Hevl  nnerP  rods 


1K  mJ  ^  =  DownDyadH i (J  Levi  nnerP  rods,  A), 


(30) 


DownDyadH  i.m  performs  the  exact  same  tasks  ( i . e.  circular  convolution,  downsampling  by  a  factor  of  2)  as 
DownDyadLo.m,  however,  it  manipulates  the  ..Iter  h  to  formulate  the  ..Iter  g  according  to  gk  =  (~l)khit  k. 
T  herefore,  just  as  was  thecase  with  thescaling  functions,  the  inner  products  can  befound  with  the  derivative 
of  the  wavelets  to  extend  the  multiresolution  down  J  levels  to  the  j  -  j  =  o  level  of  resolution,  where  the 
dyadic  decimation  will  have  left  only  1  inner  product  to  be  computed  for  the  set  WavJ_J  LevinnerProds 
[22], 


Finally,  as  we  have  seen  that  a  multiresolution  can  be  carried  out  iteratively  using  the  derivative  of  both 


a  1-D  scaling  function  and  a  1-D  wavelet,  we  recall  that  we  must  compute  the  inner  products  between  the 


derivatives  of  our  two-dimensional  wavelets  and  the  grayscale  intensity  function,  ( I 


d0\ 


and  <  I 


dd 


dx\  /  \  ’  dx2 

to  solve  Bernard’s  equations  S^.  T  he  process  we  use  to  generate  these  inner  products  relies  on  the  same 
concept  of  separability  that  was  applied  to  create  the  two-dimensional  wavelets  9n  (for  n=  l,..,4)  from  the 
tensor  products  of  one-dimensional  scaling  functions  and  wavelets,  cp  and  ip  [3],  [7],  [28], 


For  example,  let's  say  we  seek  to  ..nd  the  inner  products  ( I, 


(*1)^*2 (*2).  we  see  that 


ae). 


\ 


Ml 

dxi 


and  (  I 


deV 

’  0x2 


1.  Since d)k (xi, x2)  = 


QQ 1 

—  =  and  — ^  =  ipjkl(xi)<p°jk2{x2). 


Since  the  two- 


-jfcl  - -  T  JKl  VTJt2  V—  -  Qx 

dimensional  wavelets  are  separable  in  the  x\  and  :c2  directions,  so  are  the  2-D  inner  products.  In  other  words, 
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/  dO ^  \ 

- —)  =\l(xi),ip0jkl(x1))\l(x2),(l>jk2(x2)).  The..rst  factor  in  this  tensor  product,  VOzihV’jfci  (zi)), 

is  de..ned  so  that  we  are  taking  grayscale  intensity  values  pixel  row  by  pixel  row  (i.e.  /(£ i)  =  I(x i,x2) 
where  x2  is  ..xed)  to  generate  a  1-D  inner  product.  In  the  second  factor,  \/(x2),  <^fe2(x2)),  intensity 
values  are  being  taken  pixel  column  by  pixel  column  (i.e.  /( X2)  =  I(xi,x2)  w\iere  xi  is  ..xed).  Like¬ 


wise,  we  may  write 


dX2 


as  a  tensor  product  of  one-dimensional  inner  products, 


r 

dX2 


\l{x\), ipjkl(xi))  \/(x2),  (jJjfo (x2)).  Now,  we  see  that  we  may  apply  the  previous  analysis  for  comput¬ 
ing  the  1-D  inner  products  with  either  <j>,  ip  or  their  derivatives  at  various  levels  of  resolution  j  using  the 
cascade  algorithm.  This  provides  ail  of  the  necessary  components  for  the  construction  of  Sj k  [22],  [3],  [7], 
[28], 

As  the  ..nal  implementation  of  Bernard's  algorithm  currently  remains  an  un..nished  product,  we  turn 
to  the  completely  dinerent  approach  to  solving  for  optic  low  ..rst  proposed  by  Horn  and  Schunck  [14], 


7  Calculus  of  Variations  Approach 

7.1  An  Additional  Constraint  Equation 

Optic  low  is  no  dinerent  from  any  other  question  appearing  in  applied  mathematics  in  that  vastly 
dinerent  methods  have  been  proposed  to  ..nd  a  solution.  W  hile  Bernard  relied  on  a  multiscale  approach  to 
treat  the  problem  as  a  well-posed  matrix  equation  [3],  [4],  [5],  [6],  [7],  Horn  and  Schunck  instead  turn  to  the 
addition  of  an  additional  constraint  upon  the  motion  of  grayscale  intensity  [14],  [3], 

The  approach  is  founded  upon  the  premise  that  sharp  edges  occurring  in  an  image  sequence  cause 
competition  between  a  large  number  of  vectors  representing  the  optic  low.  In  other  words,  if  two  distinct 
objects  which  neighbor  each  other  are  moving  with  dinerent  velocities,  it  becomes  very  di<{  cult  to  determine 
the  motion  of  grayscale  intensity  along  a  sharp  boundary  between  them  [14],  T  his  concept  is  illustrated  in 
F  igure  13. 


IMAGE  1 

OPTIC  FLOW 

IMAGE  2 

■ 

1 1 1 9  1  M 

HJ! 

TTT?IU 

1 

Figure  13:  Sharp  Edges  Cause  Competition 
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Here  we  have  two  images  from  a  sequence  in  which  a  white  rectangle  borders  a  black  rectangle  with  a  gray 
background.  The  white  rectangle  moves  up  while  the  black  one  moves  down.  Describing  the  optic  fow 
poses  no  problem  except  at  the  sharp  edge  which  de..nes  the  boundary  between  the  two  objects.  Along 
this  edge,  it  becomes  unclear  whether  the  vector  ..eld  should  indicate  an  upward  or  a  downward  motion  of 
grayscale  intensity. 

According  to  Horn  and  Schunck,  this  di <t  culty  may  be  overcome  by  imposing  a  smoothness  constraint 
upon  the  optic  fow  ..eld.  In  essence,  the  constraint  carries  out  a  form  of  averaging  to  allow  for  the  ..nding 
of  one  single  "averaged"  vector  to  represent  the  optic  fow.  The  expense  of  such  a  technique  is  the  loss  of 
sharpness  [14], 


This  smoothness  constraint  can  be  expressed  as  the  minimization  of  the  square  of  the  magnitude  of  the 
gradient  of  the  optic  fow  velocity.  To  ease  the  burden  of  notation  for  the  following  discussion,  we  adopt 
x  and  y  as  the  Cartesian  coordinates  rather  than  x\  and  *2-  Furthermore,  we  refer  to  the  optic  fow  as 


having  components  u  and  v  in  v  =  (u,  v)  whereas  we  had  previously  used  v\  and  v2.  If  we  de..ne  the 


du  du 


gradients  Vtt  =  ( ux,uv )  =  j  and  Vv  =  (vx,vy)  =  —  J,  then  the  measure  of  the  departure 

from  smoothness  in  the  velocity  fow  may  be  expressed  as 


dv  dv 


4  =  !lvw||  +  ||v«||  =  v  (ux)  +  (zty) 


(VX)2  +  ( Vyf 


=  <+  Uy  +  vz  +  vr 


Furthermore,  we  incorporate  the  intensity  constraint  equation  (i.e.  the  governing  equation  for  optic  fow) 
and  represent  the  error  due  to  the  grayscale  intensity  changing  over  local  space  and  time  as 


£l  =  IXU+  IyV+  It 


.  r  di  _  di  ,  _  di 
where  7,  =  aid 

The  goal  then  becomes  to  minimize  the  total  error 


s2  =  J  J  [e2  +  <5e|)  dxdy.  (31) 

T  his  minimization  of  Equation  (31)  is  achieved  by  ..nding  the  appropriate  values  for  the  optic  low  velocity 
v  =  (u,v).  Thus  we  are  presented  with  a  scheme  for  satisfying  the  constraint  equations  by  solving  for  the 
vector  ..eld  we  seek.  We  go  about  such  a  task  using  a  method  called  Calculus  of  Variations  [14],  [18],  [13, 
Ch.  2], 
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7.2  Calculus  of  Variations 

One  of  the  key  concepts  in  Calculus  of  Variations  is  what  is  known  as  a  functional.  Consider  a  set  7  of 
functions  which  all  satisfy  certain  conditions.  A  functional  can  bethought  of  as  no  more  than  a  quantity 
which  assumes  a  speci..c  value  corresponding  to  each  function  in  7.  In  essence,  a  functional  is  a  function  of 
functions.  For  example,  the  de..nite  integral  1  =  f^f(x)dx  is  a  functional  since  its  value  is  determined  by 
the  function  /  [13,  pg.  131]. 

The  most  fundamental  problem  in  Calculus  of  Variations  is  to  ..nd  the  necessary  conditions  on  u  so  that 
it  minimizes  the  functional 


I[u) 


F  ( u ,  ux,uy)  dxdy , 


(32) 


where  F  can  be  differentiated  at  least  twice  in  both  x  and  y.  I  n  order  to  ..nd  the  necessary  conditions  on  u, 
we  consider  the  functional  J[e]  de..ned  as  I[u  +  ey\  where  u  is  the  minimizer  and  y  is  any  smooth  function 
in  the  region  -  which  vanishes  on  the  boundary  of  -  ,  denoted  by  d-  .  In  order  for  u  to  minimize  Equation 
(32),  e  =  0  must  minimize  J.  Therefore,  we  differentiate  J  with  respect  to  c  and  set  it  equal  to  0  when 
e  =  0.  However, 


J[e\ 


,  ux  +  e?7 


X 


uv  +  erly)  dxdy, 


(33) 


which  yields 


J°[0]  = 


II  (vFu  +  VxFu 


+  VyKy)  dxdy  =  0, 


(34) 


where  Fu  =  —,  FUx 

ou 

We  now  note  that 


dF 

=  — ,  and  Fu 

aux 


dF 

dUy 


VxFux  —  v(Fux)x  +  ( vFUx)xi  Vx  FUx  —  V(FUy)y  +  ( f]FUy)y 


If  we  apply  the  de..nition  of  divergence,  div(m,n) 


dm 

dx 


— ,  we  see  that 

dy 


A]  =  1 1  V  (FU ( FUx)x  +  (FUy)y)  dxdy  +  ^  jT  div  (y\FUx,  FUy ))  dxdy  =  0. 


(35) 


Next  we  apply  the  Divergence  Theorem  which  states:  J  j  div(\N  )dxdy  =  I  (W  •  n)  dt  where  n  is 
the  outward  normal  vector  to  d-  .  This  allows  us  to  replace  the  second  integral  in  Equation  (35)  with 
I  ((v\Fux, FUy))  ■  n)  dt.  Since  y  was  de..ned  to  satisfy  zero  boundary  conditions  on  d-  ,  this  integral 
vanishes,  \herefore,  we  may  write 
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J°[°]  =  J  j  V  (pu  ~  ( FUJX  +  ( FUy)v )  dxdy  =  0,  ( 36) 

which  applies  for  any  smooth  function  77  in  -  .  In  order  for  this  to  be  the  case,  the  integrand  must  itself  be 
equal  to  zero.  T  his  provides  the  necessary  condition  on  F  and  u  for  there  to  be  a  minimizer: 

Fu  -  (FUx)x  -  (FUy)v  =  0-  (37) 

Equation  (37)  is  referred  to  asthe  E uler-Lagrange  equation  for  Equation  (32)  [18],  [13,  Ch.  2], 

Let  us  now  examine  an  example  problem,  namely  Laplace's  Equation.  Consider  the  Dirichlet  integral 

./[u]  =  J  J  (u2  +  Uy)  dxdy.  (38) 

Here  we  see  that  F(u,ux,uy)  =  u2x  -fu*.  We  now  apply  Equation  (37), 

Fu  ~  (FUx)x  —  ( FUy)y  =  0 

0  (2x^3;)^  [2'Uy)  y  -  0 

“L  yy )  -  0 

So  that  we  obtain 


V2zt  =  0, 


(39) 


where  \/2u  =  uxx  +  uvv  =  +  — rr  is  the  Laplacian  of  u  [18],  [13,  C  h.  2], 

ox1  ayz 

We  now  seek  to  show  that  the  same  argument  may  be  applied  to  a  functional  J  which  operates  on  a 
vector  u  =  (zt,  v ).  I  n  this  case,  we  seek  to  minimize 


J[u ,  v] 


F(u,  v ,  ux,  uy ,  vx,  vy)dxdy. 


Applying  the  principle  of  the  Euler-Lagrange  equation  to  two  dimensions  we  see  that 


(40) 


Fu  ~  (FUx)x  -  (FUy)y  =  0,  Fv-  (FVx)x- (FVy)y  =  0.  (41) 

We  now  have  the  basic  tools  needed  to  apply  Calculus  of  Variations  to  the  two-dimensional  context  of 
optic  fow  stated  as  a  minimization  problem  in  Equation  (31)  [18],  [13,  Ch.  2], 
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7.3  Application  of  Calculus  of  Variations  to  Optic  Flow 

Recalling  our  de..nitions  of  £7  and  s2,  we  rewrite  Equation  (31)  as  a  functional  with  arguments  u  and 
v  representing  the  components  of  the  optic  tow,  v  =  (u,v). 

J[u,v]  =  j  J  (It  +  Ixu+  Iyv)2  +  5(u2  +  u2  +  v2  +  v2)dxdy.  (42) 

Since  theterm  S(ul  +  u2  +  vl  +v2)  represents  thedeparture  from  smoothness  in  the  velocity  tow,  it  may 
be  viewed  as  making  the  problem  less  "rigid"  mathematically  in  that  it  allows  for  a  certain  amount  (i.e.  an 
amount  proportional  to  the  regularization  parameter  5)  of  elastic  stretching  in  the  ..eld  which  de..nes  optic 
fow  [1], 

We  now  apply  Equation  (41)  and  the  result  obtained  from  Laplace's  Equation, 

Fu  ~  (Fux)x  +  (FUy)i i  —  d,  Fv  —  {FVJ)X  +  (FVy)y  =  0 
2  Ix(It  +  Ixu  +  Iyv)  —  <5(2V2u)  =  0,  +  Ixu  +  Iyv)  —  S(2\72v)  =  0 


which  yields 


Ix(It  +  Ixu  +  Iyv)  =  SV2u,  Iy(It  + Ixu  +  Iyv)  =  SV2  V.  (43) 

T  hese  are  the  same  equations  which  Horn  and  Schunck  attempt  to  solve  for  the  optic  fow  v  =  (u,v)  [18], 
[14],  Notice  that  there  are  two  equations  with  two  unknowns,  making  the  problem  well-posed. 

At  this  point,  our  method  diverges  from  that  of  Horn  and  Schunck.  We  seek  to  apply  the  Galerkin 
method  in  order  to  ..nd  a  numerical  solution  to  the  system  of  partial  di nerential  equations  in  (43)  whereas 
an  approach  based  on  ..nite  dinerences  is  proposed  in  [14], 

7.4  Insertion  of  the  Galerkin  Method 

Our  departure  from  the  ..nite  dinerence  scheme  presented  by  Horn  and  Schunck  bears  no  relation  to  a 
lack  of  trust  in  their  method.  On  the  contrary,  the  theoretical  development  behind  ..nite  dinerences  shows 
strong  promise  for  application  to  such  a  problem  as  optic  tow  [14],  [10],  [25],  [15], 

Rather,  we  choose  the  Galerkin  method  based  on  the  fact  that  it  seeks  solutions  of  partial  dinerential 
equations  in  terms  of  a  countable  basis  <pt  [18],  Extensive  work  has  been  done  to  prove  the  enectiveness 
of  the  Galerkin  or  "spectral"  method  in  solving  PDEs  using  traditional  bases  such  as  Fourier  (sines  and 
cosines)  or  polynomials  (Chebyshev,  Legendre,  or  Hermite)  [11],  [21],  [12],  [25],  [9],  [18],  Thus  the  appeal 
of  applying  the  Galerkin  method  to  optic  tow,  aside  from  the  fact  that  it  has  never  been  done  before,  is 
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that  it  lends  itself  to  the  employment  of  a  wavelet  basis.  Although  the  strategy  for  making  use  of  wavelets 
in  this  context  diners  greatly  from  B  ernard's  approach,  there  exists  the  possibility  that  they  will  open  up  a 
new  door  in  the  overall  quest  to  solve  the  problem  of  optic  fow.  For  now,  we  describe  the  method  in  terms 
of  a  generic  countable  basis  fc. 

Going  back  to  Equation  (36),  we  choose  a  special  class  of  functions  rj,  namely  4>t,  with  the  end  goal 
being  to  seek  solutions  of  the  form 

11  ii 

u(x,y,t)  =  ^2^2ai,  ,  v{x:y,t)  =  <44) 

i=i>=i  J=lt=l 

Here,  we  see  that  the  components  of  optic  fow  can  be  written  as  an  exact  sum  of  time  dependent  coep  cients 
multiplying  a  two-dimensional  spatially  dependent  basis.  In  practice,  we  are  only  able  to  compute  a  ..nite 
number  of  terms  in  Equation  (44).  If  we  use  enough  terms,  then  the  equations  in  (44)  can  be  approximated 

N  N  N  N 

u(x,  y,  t)  «  uN  (x,  V,t)  =  a %  j  (t)^(a;)^  (y),  v(x,  y,  t)  «  vN(x,  y,  t)  =  ^  ^  bit  j  (i)<4(4<4  (y). 

j  =  1  2  =  1  j  =  l  i=l 

(45) 


Now,  applying  Equations  (36)  and  (43),  we  see 


II  (K  -  (FUx)x  -  ( FUy)y )  |„= 


UN,  V  =  VN 


dxdy  =  0,  *=1,2,  ...,7V  and  j  =  l,2,...,N 


1 1  (4  (4  +  hu  +  Iyv)  -  SV2u)  \u=un,  v=vn  dxdy  S  0, 


i  =  1,  2, ...,  N  and  j  =  1,  2, ...,  N 


(46) 


and 


/  I.  44)  (Fv  (FVx)x  ( FVy)y  )  |u  = 


u  N,  v=  VN 


dxdy  —  0,  7  =  1,2,  ...,7V  and  j  =  l,2,...,N 


1 1  (4  (4  +  4«  +  V)  -  *v2«)  L= 


uNi  V  =  VN 


dxdy  =  0,  i  =  1,  2, ...,  N  and  j  =  1,  2, ...,  N. 


(47) 


Equations  (46)  and  (47)  represent  2N2  equations  with  2N2  unknowns.  Thus,  we  have  utilized  the 
smoothness  constraint  equation,  calculus  of  variations,  and  the  Galerkin  method  to  create  a  solvable, 
well-posed  problem  [18],  [14],  The  di<t  cuity  now  lies  in  ..nding  a  feasible  scheme  for  implementation  in 
M  athematicaw ,  our  scienti..c  computing  program  of  choice  for  this  method. 


42 

7.5  Issues  Of  Implementation 

In  both  Equations  (46)  and  (47),  we  observe  the  presence  of  the  term  It,  the  time  derivative  of  the 
grayscale  intensity  function.  The  de.nition  of  this  derivative  states  that 

r  dl(x,y,t)  I(x,y,t  +  h)~  I(x,y,t) 

at  h\  o  h 

A  problem  arises  as  one  tries  to  implement  the  Galerkin  approach  to  solve  for  optic  tow  in  that  the  value 
of  h  cannot  be  taken  to  the  limit  of  0.  The  time  step  is  ..xed  by  the  interval  between  discrete  frames  in 
the  image  sequence.  Thus  we  must  choose  an  approximation  to  For  our  purposes,  we  have  very  simply 
chosen  h  =  1  so  that  It  =  I(x,y,t  +  1)  -  I(x,y,t).  Although  this  is  dinerent  from  the  approximation  given 
in  [14],  no  evidence  has  presented  itself  so  far  to  show  that  our  approach  will  speci . .cal ly  cause  instability. 

Another  issue  presents  itself  with  relation  to  the  spatial  derivatives  of  the  intensity  function,  Ix  and 
Iy.  In  order  to  take  these  derivatives  directly,  we  must  have  a  function  I  which  is  continuous  in  the  x 

and  y  directions.  This  is  accomplished  by  interpolating  the  discrete  grid  of  pixel  intensity  values  obtained 

(R) 

by  importing  a  grayscale  image  into  M  athematica  .  The  command  which  performs  this  task  is  called 
Listl  nterpolati  on. 

The  other  spatial  derivatives  being  taken,  namely  the  Laplacians  V2w  and  X72v,  do  not  pose  such  an 

obstacle  because  the  components  of  the  velocity  can  be  represented  in  terms  of  dinerenti able  basis  functions, 

ii  ii 

u{x,y,t)  =  EE  v(x,  y,  t)  =  EE  K 

j  =  H=  l  j  =  H  =  l 

D inerenti abi lity  of  the  cosine  function  is  one  reason  for  our  choice  of  the  Fourier  basis  to  serve  as  <j>. 
T  he  other  reasoning  comes  from  the  fact  that  there  already  exists  a  solid  foundation  for  use  of  this  basis  in 
the  Galerkin  method  [11],  [21],  [12],  [25],  [9],  [18],  There  exists  a  de..nite  problem,  however,  in  applying 
periodic  bases  of  in..nite  support,  such  as  cosines,  to  functions  which  are  of  ..nite  domain  and  are  most  likely 
non-periodic,  such  as  the  optic  fow  generated  by  a  grayscale  image  sequence.  Traditionally,  to  attain  a 
good  quality  approximation  of  such  a  function  using  the  basis  {<fin(x)  =  cos(n7ra)}n,  a  very  large  number 

of  modes  are  required.  This  means  that  the  "cut  on"  number  N  will  have  to  be  very  large  for  u  «  uN 

(r) 

and  v  «  vN.  Thus  our  computation  of  the  optic  fow  vectors  which  requires  M  athematica  to  solve  2 N2 
integral  equations  for  2 N2  unknowns  will  take  an  extremely  long  amount  of  time  and  may  even  exhaust  the 
memory  for  a  standard  computer  [12],  [26], 

As  was  mentioned  earlier,  the  Galerkin  method  can  be  executed  using  a  varidty  of  dinerent  choices 
for  basis  functions  cp.  T he  ditf  culty  of  approximating  the  optic  fow,  as  described  in  the  preceding  para¬ 
graph,  would  be  be  better  approached  using  a  basis  whose  functions  were  indexed  over  both  space  and 
scale,  such  as  wavelets  [26],  [16],  The  only  problem  with  using  such  a  basis  is  that  the  approxima- 

N  N  N 

tions  to  the  optic  4 ow  then  become  uN(x,y,t)  =  k2(t)4>j,  ki(x)<t>j,k2(y)  and  vN(x,y,t)  = 

j  =  1  ki  =  1  k2  =  1 
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N  N  N 

^2^2  ^2  bj ,  ki,  ki ki(x)<t>j,  k2(y)  which  when  substituted  into  Equations  (46)  and  (47)  result  in  the 

j=l ki=l fc2=l 

solving  of  2  N3  equations  for  2  N3  unknowns.  Consequently,  the  advantage  of  requiring  less  terms  for  a  good 
approximation  may  not  be  gained  after  all. 

As  is  the  case  at  the  completion  of  every  research  project,  a  number  of  issues  still  remain  unresolved. 
T  herefore,  let  us  now  turn  to  the  results  obtained  using  the  protocol  for  implementation  which  is  already  in 
place. 


8  Results 

8.1  Wavelet  Based  M  ultiscale  Approach 

Results  have  not  yet  been  attained  using  this  approach.  The  remaining  unresolved  issues  of  implemen¬ 
tation  are  extremely  close  to  being  overcome. 

8.2  Calculus  of  Variations  A pproach 

The  ..rst  attempt  at  ..ndingan  optic  fow  vector  ..eld  wasconducted  with  a  very  simpli..ed  approximation 
to  motion  occurring  between  two  images.  Since  a  grayscale  image  is  treated  as  a  surface  by  a  computer,  we 
simply  used  known  functions  to  create  surfaces  that  would  represent  two  pseudo-images.  Our  ..rst  image  in 
the  sequence  was  represented  by  the  function  el(x,  y)  =  cos(27rx)  sin(37ry)  over  the  domain  [0,1]  x  [0,1].  In 
order  to  simulate  motion  we  add  to  el  a  bump,  and  so  create  the  second  image  represented  by  the  function 
e2(x,y)  =  el (x,y)  +ga,b{x,y),  where 

gab(x,y)~ s.  1  cos^  '  ra,b(x,y))  if  ra,b(x,y)  =  s/(x  -  a)2  +  {y  -  b)2  <  - 
[  0  otherwise 

and  we  set  a  =  0.45  and  b  =  0.53. 

Figure  14  illustrates  the  ..rst  and  second  pseudo-images  and  then  the  bump  which  simulates  the  motion 
of  grayscale  intensity  between  the  two. 


Pseudo-Image  1  Bump  Representing  Pseudo-Image  2 

Motion 


Figure  14:  Simple  Case  Pseudo-Images  and  Bump 
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(R) 

T  he  corresponding  optic  fow  vector  ..elds  were  solved  for  with  the  M  athematicaw  notebook  HSl.nb  using 
both  3  and  11  modes  for  the  Fourier  cosine  basis.  Figure  15  displays  the  results. 

Optic  Flow  with  3  Modes  Optic  Flow  with  1 1  Modes 


F  igure  15:  Optic  F  low  B etween  Simple  Case  Pseudo-Images 


For  both  cases,  there  are  non-trivial  vectors  occurring  only  in  the  region  surrounding  x  =  .45  and  y  =  .53. 

(R) 

T  he  rest  of  the  arrows  are  represented  by  M  athematica  as  having  no  tails  and  so  they  can  be  assumed  to 
be  of  value  zero.  T  his  result  is  as  it  should  be  since  the  motion,  simulated  by  the  bump,  only  occurred  in  a 
region  centered  at  (.45,  .53).  Looking  closer  at  the  vector  ..elds,  wesee  that  the  motion  appears  to  betaking 
place  over  a  larger  area  when  we  used  3  modes  to  solve  for  it  as  opposed  to  using  11.  If  we  glance  back 
to  Figure  14,  we  see  that  the  relative  size  of  the  bump  in  the  domain  [0, 1]  x  [0, 1]  more  closely  resembles 
the  size  of  the  region  experiencing  motion  in  the  optic  fow  ..eld  with  a  solution  based  on  11  modes.  T  his 
provides  evidence  to  support  our  conjecture  that  the  approximation  to  the  optic  fow  is  more  accurate  if  we 
use  more  basis  functions.  The  only  problem  is  that  the  solution  took  signi ..cantly  more  time  ..nd  when 
M  athematica®  had  to  solve  the  2(ll)2  =  242  equations  instead  of  the  2(3)2  =  18  equations  (approximately 
30  minutes  as  opposed  to  2  minutes). 

We  see  that  the  basic  premise  of  optic  fow  has  been  put  into  place  for  a  sim pi i . . ed  approximation  to 
an  image  sequence.  T  he  next  step  is  to  ..nd  solutions  to  the  optic  fow  bdtween  two  grayscale  images  in  an 
actual  sequence.  T  his  was  accomplished  using  a  slight  variation  from  our  method  described  earlier. 
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Previously,  we  used  the  command  L i stl nter pol ati on  to  create  functions  which  represented  the  images  so 

that  we  could  take  the  partial  derivatives  lx  and  Iv.  The  presence  of  these  interpolated  functions  is  one  of 

(?) 

the  leading  factors  which  cause  M  athematica  to  take  a  long  time  to  solve  the  integral  equations  in  (46) 
and  (47).  T  herefore,  we  attempt  to  circumvent  this  problem  by  using  a  Fourier  cosine  series  to  approximate 
the  images  according  to 


N  N 

I(x,y)  ~  IN(x,y)  =  E  E  Cm.n<f>m(2:)0n(2/), 

n=ln=l 


where  the  coe<t  cients  are  computed  by  cm  „  =  T7- 

T  hus,  despite  ail  of  the  partial  derivatives  and  basis  functions  contained  in  (46)  and  (47),  the  integrands  of 
these  equations  become  nothing  more  than  a  multiplication  between  sines  and  cosines.  Solving  a  system 
such  as  this  is  much  less  computationally  intensive  since  the  integrals  can  be  computed  exactly  without  the 
need  to  consider  the  ejects  of  interpolation.  The  drawback  to  this  technique  is  that  accuracy  is  lost  in 
performing  a  cosine  approximation  of  an  image.  Sharp  transitions  in  grayscale  intensity  are  di<t  cult  to 
capture  using  functions  which  are  smooth  and  of  in.mite  support.  The  bene..t,  however,  of  this  approach  is 
that  we  are  able  to  use  more  frequency  modes  in  the  approximation  while  maintaining  the  time  required  to 
..nd  a  solution  to  the  optic  fow. 

Figure  16  illustrates  the  fact  that  the  general  trend  of  motion  can  be  captured  with  this  method  using 

(?) 

16  modes  of  frequency.  The  optic  fow  is  found  using  the  M  athematica  notebook  H SG aler ki n. n b. 


F  igure  16:  Optic  F  low  U  sing  the  Galerkin  M  ethod 


9  Conclusions 
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We  have  explored  some  of  the  fundamental  theory  which  lies  behind  two  dinerent  methods  of  solving  for 
optic  fow.  Additionally,  a  solid  path  has  been  laid  down  for  the  process  of  implementing  each  technique. 
M  any  riddles  and  obstacles  that  existed  before  have  now  been  resolved  or  at  least  have  a  starting  point  from 
which  to  approach  them. 

It  remains  inconclusive  as  to  which  is  a  better  method,  the  wavelet  based  multiscale  approach  proposed 
by  Bernard  or  the  calculus  of  variations  scheme  de/eloped  by  Horn  and  Schunck.  On  one  hand,  Bernard's 
technique  appears  to  be  more  robust  in  its  ability  to  address  the  issues  of  time  aliasing  and  aperture.  T  he 
multiresolutional  nature  of  the  wavelet  basis  seems  to  make  it  the  perfect  tool  for  dealing  with  such  a 
problem.  A  weakness  arises  in  Bernard's  argument,  however.  He  makes  the  assertion  that  by  carrying 
out  the  multiresolutional  decomposition  all  the  way  to  the  coarsest  I  e/el,  the  solution  to  the  system  of 
matrix  equations  will  converge  to  the  actual  optic  fow.  The  problem  stems  from  his  key  assumption  that 
the  optic  fow  is  constant  over  the  support  of  the  wavelets  [3],  [4],  [5],  [6],  [7],  At  the  coarser  levels  of 
resolution,  the  wavelets  have  support  which  span  a  large  portion  of  the  image.  T  hat  the  optic  fow  would 
be  uniform  throughout  so  large  an  area  is  not  a  realistic  statement.  Yet  without  this  assumption  in  place, 
the  matrix  equations  cannot  be  constructed.  Therefore,  from  a  sense  of  practicality,  we  must  end  the 
multiresolutional  decomposition  at  the  level  where  the  optic  fow  begins  to  violate  this  assumption.  T  he 
question  as  to  whether  or  not  the  solution  of  this  abridged  system  will  still  converge  to  the  true  optic  fow 
is  left  unanswered  by  Bernard.  We  see  that  here  he  has  left  an  enormous  distance  between  the  theory  and 
the  practical  implementation  of  his  method. 

On  theother  hand,  the  approach  of  Horn  and  Schunck  appears  to  make  more  sense  from  a  standpoint  of 
physics.  The  key  assumption  made  by  Bernard  about  the  constancy  of  optic fow  over  various  local  scales  is 
replaced  with  the  much  more  plausible  notion  that  the  optic  fow  velocity  makes  up  a  smooth  ..eld.  Consider 
F  igure  13  again  where  two  neighboring  objects  are  undergoing  dinerent  motion.  T  he  smoothness  constraint 
of  Horn  and  Schunck  forces  the  optic  fow  to  take  on  an  "averaged"  value  at  the  boundary  between  the 
two  objects,  thus  losing  the  quality  of  sharpness  [14],  Such  a  loss  of  information  can  be  considered  much 
less  si gni . .cant,  however,  than  that  experienced  by  having  to  classify  the  motion  as  one  constant  value  for 
the  area  taken  up  by  the  entirety  of  both  objects  as  would  happen  under  Bernard's  method.  Problems 
are  encountered  with  the  calculus  of  variations  algorithm,  however,  when  the  image  sequence  is  complicated 
by  the  appearance  of  numerous  small  objects  travelling  with  distinct  velocities  either  in  close  proximity  or 
overlapping  each  other.  In  a  case  such  as  this,  the  loss  of  sharpness  experienced  at  the  borders  of  the 
objects  is  much  more  substantial  since  the  objects  may  not  comprise  su<t  cient  pixel  area  to  preserve  their 
true  motion  from  the  averaging  requirement  imposed  by  the  smoothness  constraint.  T  he  problem  with  this 
fact  is  that  almost  any  video  sequence  could  realistically  contain  this  kind  of  complexity.  T  hus,  we  see  that 
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this  method  may  only  have  practicality  for  simple  image  sequences. 

Of  course,  the  development  of  each  mdthod  has  not  reached  a  ..nal  end.  Numerous  opportunities 
exist  to  adjust  each  model  to  better  suit  the  problem  of  optic  fow.  For  instance,  Bernard's  research  has 
begun  to  explore  the  concept  of  "image  warping".  This  entails  estimating  the  dynamic  of  motion  as  a 
uniform  translation  across  the  image  (i.e.  all  objects  travel  with  the  same  velocity)  and  then  solving  for  the 
deviation  from  this  warped  image  at  various  scales.  The  assertion  is  that  if  enough  scales  are  used  in  the 
multiresolution,  the  sum  of  the  large  motion  estimates  with  the  residual  deviations  will  be  suit  ciently  close 
to  the  true  motion  [7],  As  was  mentioned  earlier,  the  use  of  the  Galerkin  method  to  solve  the  calculus  of 
variations  problem  developed  by  Horn  and  Schunck  might  be  improved  if  a  more  high-powered  basis  was 
used  such  as  one  incorporating  wavelets.  It  is  a  distinct  possibility  that  there  exists  some  trick  which  relies 
on  the  orthonormality  of  wavelets  to  make  the  solving  of  the  integral  equations  replaceable  by  a  much  less 
computationally  intensive  task  [28], 

In  addition,  some  of  the  fundamental  issues  relating  to  optic  fow  are  left  unresolved.  For  instance,  let 
us  say  that  one  was  attempting  to  track  the  real  motion  of  objects  by  ..Iming  them  and  then  solving  for 
the  corresponding  optic  fow.  How  does  one  describe  the  dinerence  between  the  actual  and  the  projected 
motion?  The  method  for  overcoming  this  obstacle,  as  mentioned  earlier,  is  presented  in  [23]  as  an  issue  of 
stochastic  error.  Si  moncell  i  states  that  it  is  possible  to  use  an  uncertainty  model  to  express  the  departure 
of  the  optic  fow  ..eld  from  the  real  world  scene  of  motion.  It  may  well  prove  that  this  approach  does  a 
reasonably  good  job  of  addressing  the  problem.  A  highly  accurate  and  adaptive  stochastic  model  would 
have  to  be  developed  ..rst,  however,  to  apply  this  idea  of  motion  sensing  based  on  optic  fow  to  a  context 
relevant  to  the  military  such  as  target  tracking. 

Another  si gn i . . cant  yet  unresolved  issue  is  that  the  basic  physics  behind  the  optic  fow  model  might 
not  hold.  It  is  possible  that  the  theory  behind  the  governing  equation  (i.e.  that  grayscale  intensity  is 
conserved)  is  not  valid.  Bernard  makes  an  attempt  to  address  this  concern  by  proposing  a  Lambertian 
surface  aspect  model  to  take  into  account  illumination  changes  [3],  Although  this  idea  makes  intuitive  sense 
as  an  improvement  on  the  way  to  solve  for  optic  X ow,  it  is  not  directly  stated  in  his  publications  that  better 
results  were  attained  through  the  use  of  it. 

Overall,  it  may  be  stated  that  optic  fow  has  a  long  way  to  go  before  it  may  stand  alongside  some  of 
the  problems  in  applied  mathematics  which  are  regarded  as  high  powered  and  crucial  for  solving.  The 
approaches  taken  to  resolve  the  issues  have  only  been  explored  to  an  extent  where  optic  X ow  may  be  applied 
in  a  limited  sense.  Questions  of  theory  are  not  met  with  answers  of  implementation  across  a  wide  range  of 
issues.  This  project  has  closed  some  of  the  gaps  in  the  pursuit  of  providing  more  potential  for  application 
of  optic  fow.  However,  much  work  still  remains  for  future  studies  to  build  higher  upon  the  foundations.  It 
will  be  interesting  to  observe  how  far  research  is  able  to  advance  the  concept  and  what  it  will  be  applied  to. 


It  is  the  hope  of  the  author  that  some  day  enough  will  be  known  to  harness  optic  fow  for  use  in  the  domain 
of  acoustic  imagery  supplied  by  sonar  systems. 
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11  Appendix 


MAT  LAB®  m-..les 


Haarl.m 

%Used  to  compute  Level-1  Haar  Transform 
fun cti  on  [al  evel  l,d  I evel  1]=  H  aarl(  f ) ; 

%lnput  f  is  the  sample  function  to  be  decomposed 
N  =length(f); 

%lndex  j  is  the  scale  or  level 

j=i; 

%Creation  of  matrix  W  to  store  1-level  Haar  wavelets 

W=zeros(N/2,N); 

for  m=  1:  N  /  2 

W(m,m*2/'j-(2/'j-l):m*2/'j)  =  [ones(l,2/'(j-l))/(2/'(j/  2)),-l*ones(l,2/v(j-l))/(2/v(j/  2))]; 
end 

%Creation  of  matrix  V  to  store  1-level  Haar  scaling  functions 
V  =zeros(N/  2,N ); 
for  m=  1:  N  /  2 

V (m,m*2/vj-(2/vj-l):m*2''j)  =  ones(l,2'vj)/  (2/v(j/  2)); 
end 

%F  inding  vectors  alevell=al,a2 . aN/  2  and  dlevell=dl,d2,...,dN/  2 

for  m=  1: N /  2 

a(m)=f*V(m,:)’; 

d(m)=f*W(m,:)'; 

end 

alevell=a; 


dl  evel  l=d ; 


HaarlO.m 

%Used  to  compute  Level-10  Haar  Transform  &  plot  Averaged  Signals 
function  A  =  HaarlO(f); 

%l  nput  f  is  the  sample  function  to  be  decomposed 
N  =length(f); 

%Creation  of  matrix  W  to  store  Haar  wavelets 
W  =zeros(N/  2,N  ,10); 
for  j  =  1:10 

for  m  =  l:N/(2/vj) 

W(m,m*2/'j-(2/'j-l):m*2/'j,j)  =  [ones(l,2/'(j-l))/  (2~(j/  2)),-l*ones(  l,2^(j-l))/  (2~(j/ 2))  ]; 
end 
end 

%Creation  of  matrix  V  to  store  Haar  scaling  functions 
V  =zeros(N/  2,N  ,10); 
for  j  =  1:10 

for  m  =  l:N/(2/vj) 

V(m,m*2~j-(2~j-l):m*2~jJ)  =  ones(l,2~j)/(2~(j/2)); 
end 
end 

%F  in di ng  vectors  alevelj  and  dlevelj 
for  j  =  1:10 

for  m  =  l:N/(2/vj) 

a(m,j)=f*V  (m,:,j)’; 
d(  m,j)=f*W(m,:,j)'; 
end 
end 

%F  inding  the  ..rst  10  averaged  signals  A 1  to  A10 

A=zeros(10,N); 

for  j  =  1:10 

temp=zeros(l,N ); 
for  m  =  l:N/(2/vj) 

temp=temp+a(m,j)!t!V  (m,:,j ); 
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end 

A(j,:)  =  temp; 
end 

%T emporary  code 

%atemp=a(..nd(a(:,2)~=0),2);  dtemp=d(..nd(d(:,2)~=0),2); 

%A  =  [atemp',dtemp',d(:,l) ']; 

%T emporary  code 
alOtemp=a(  1:2^0,10); 

dl0temp=d(l:2~0,10);  d9temp=d(l:2/vl,9);  d8temp=d(l:2/v2,8); 
d7temp  =  d(l:2/v3,7);  d6temp=d(l:2/v4,6);  d5temp=d(l:2/'5,5); 
d4temp  =  d(l:2/'6,4);  d3temp=d(l:2/v7,3);  d2temp=d(l:2~8,2); 
A=[al0temp',dl0temp',d9temp',d8temp',d7temp',d6temp',d5temp', 
d4temp’,d3temp',d2temp’,d(:,l)’]; 

Daub4_  10.  m 

%Used  to  compute  Level-10  Daub4  T ransform  &  plot  Averaged  and  Detail  Signals 
function  [A,D]=Daub4_  10(f); 

%l  nput  f  is  the  sample  function  to  be  decomposed 
N  =length(f); 

%Creation  of  matrix  V  to  store  Daub4  scaling  functions 
V  =zeros(N/  2,N  ,10); 

%Scal in g  N  umbers  alphal,  alpha2,  alpha3,  alpha4 

alphal=  (l+sqrt(  3) ) / (4* sqrt( 2) ) ;  al ph a2= ( 3+ sqrt( 3) )/  (4*sqrt(2)); 

alpha  3=  ( 3-  sq  r  t  ( 3) )/  (4*sqrt(  2)) ;  al  ph  a4= ( l-sqrt(  3) )/  (4*sqrt(2)); 

%C  reate  1st  level  scaling  functions,  then  from  this  generate  the  next  levels  in  the  M  RA 
V(l,l:4,l)  =  [alphal,alpha2,alpha3,alpha4]; 
for  m=  2: ( N /  2) 

V(m,:,l)  =  [V(l,N-2*(m-l)+  l:N,l),V(l,l:N-2*(m-l),l)]; 
end 

%C  reate  2nd  through  10th  level  scaling  functions  from  1st  level  scaling  functions 
for  j  =  2:10 

for  m=l:N/(2/vj)-l 


V  (m,:,j)  =  alphal*V  ( 2*m- 1,  :,j -1)  +  alp  ha2*V  (2*m,:,j-l)  + 
alpha3*V(2*m  +  l,:,j-l)  +  alpha4*V(2*m  +  2,:,j-l); 
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end 

%Account  for  wrap  around  [V_  2m+l&  V_  2m+2  will  exceed  N/(2/vj)] 
m=N/  (2~j ); 

V(m,:,j)  =  alphal*V  ( 2*m-l,:,j-l)  +  alpha2*V(2*m,:,j-l)  + ... 
alpha3*V(2*m+l-N/  (2/v(j-l)),:,j-l)  +  alpha4*V (2*m+2-N/  (2/'(j-l)),:,j-l); 
end 

%Creation  of  matrix  W  to  store  Daub4  Wavelets 
W  =zeros(N/  2,N  ,10); 

%Scal in g  N  umbers  betal,  beta2,  bdta3,  beta4 

betal=(  1- sqrt(3) )/  (4*sqrt(2));  beta2=  (sqrt ( 3) -3)/  ( 4* sqrt(  2) ) ; 

beta3=(3+sqrt(3))/  (4*sqrt(2));  beta4=(-l-sqrt(3))/  (4*sqrt(2)); 

%C  reate  1st  level  Wavelets,  then  from  this  generate  the  next  levels  in  the  M  RA 
W  ( 1, 1:4, 1)  =  [beta  I,beta2,beta3,beta4]; 
for  m=  2:( N /  2) 

W(m,:,l)=  [W  (l,N-2*(m-l)  +  l:N,l),W  (l,l:N-2*(m-l),l)]; 
end 

%C reate  2nd  through  10th  level  Waveletsfrom  1st  level  Wavelets 
for  j  =  2:10 

for  m  =  l:N/(2/vj)-l 

W  (m,:,j)=  betal*V (2*m-l,:,j-l)  +  beta2*V (2*m,:,j-l)  +  ... 
beta3*V(2*m+  l,:,j-l)  +  beta4*V(2!)!m+  2, :  ,j  - 1) ; 
end 

%Account  for  wrap  around  [W_  2m  +  l  &  W_  2m  +  2  will  exceed  N/(2/'j)] 
m=N/(2/vj); 

W(m,:,j)  =  betal*V  (2*m-l,:,j-l)  +  beta2*V(2*m,:,j-l)  +  ... 
beta3*V(2*m  +  l-N/  (2/v(j-l)),:,j-l)  +  beta4*V(2*m  +  2-N/(2/v(j-l)),:,j-l); 
end 

%F  in di ng  vectors  alevelj  and  dlevelj 
for  j  =  1:10 

for  m  =  l:N/(2/vj) 
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a(m,j)=f*V  (m,:,j)'; 
d(  m,j)=f*W(m,:,j)'; 
end 
end 

%F  inding  the  ..rst  10  averaged  signals  A 1  to  A10 
A  =zeros(10,N ); 
for  j  =  1:10 

temp=zeros(l,N ); 
for  m  =  l:N/(2/vj) 

temp=temp+a(m,j)*V  (m,:,j); 
end 

A(j,:)  =  temp; 
end 

%F inding  the  ..rst  10  detail  signals  D1  to  DIO 
D=zeros(  10,N ); 
for  j  =  1:10 

temp=zeros(l,N ); 
for  m  =  l:N/(2/vj) 

temp=temp+d(m,j)*W  (m,:,j); 
end 

D(j,:)  =  temp; 
end 

MakeDaub4_  10.  m 

%Purpose  is  simply  to  make  various  Daub4  scaling  functions  and  wavelets  on  their  own 
function  [V,W  ]=M  akeDaub4_  10; 

%Arbitrarily  set  N  =  1024  since  many  of  the  signals  we  have  been  using  are  this  length 
N  =2~10; 

%Creation  of  matrix  V  to  store  Daub4  scaling  functions 
V  =zeros(N/  2,N  ,10); 

%Scal in g  N  umbers  alphal,  alpha2,  alpha3,  alpha4 
alphal=(l+sqrt(3))/(4*sqrt(2));  al ph a2= ( 3+ sqrt( 3) )/  (4*sqrt(2)); 
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alpha3=(3-sqrt(3))/(4*sqrt(2));  al pha4= ( l-sqrt( 3) )/  (4*sqrt(2) ); 

%C  reate  1st  level  scaling  functions,  then  from  this  generate  the  next  levels  in  the  M  RA 
V(l,l:4,l)  =  [alphal,alpha2,alpha3,alpha4]; 
for  m=  2:( N /  2) 

V(m,:,l)  =  [V(l,N-2*(m-l)+  l:N,l),V(l,l:N-2*(m-l),l)]; 
end 

%C  reate  2nd  through  10th  level  scaling  functions  from  1st  level  scaling  functions 
for  j  =  2:10 

for  m  =  l:N/(2/vj)-l 

V(m,:,j)  =alphal*V(2*m-l,:,j-l)+alpha2*V(2*m,:,j-l)  +  ... 
alpha3*V(2*m  +  l,:,j-l)  +  alpha4*V(2*m  +  2,:,j-l); 
end 

%Account  for  wrap  around  [V_  2m+l&  V_  2m+2  will  exceed  N/(2/vj)] 
m=N/  (2~j ); 

V(m,:,j)  =  alphal*V  ( 2*m-l,:,j-l)  +  alpha2*V(2*m,:,j-l)  + ... 
alpha3*V(2*m+l-N/(2~(j-l)),:,j-l)+alpha4*V(2*m+2-N/(2/'(j-l)),:,j-l); 
end 

%Creation  of  matrix  W  to  store  Daub4  Wavelets 
W  =zeros(N/  2,N  ,10); 

%Scal in g  N  umbers  betal,  beta2,  bdta3,  beta4 

betal=(  1- sqrt(3) )/  (4*sqrt(2));  beta2=(sqrt(3)-3)/(4*sqrt(2)); 

beta3= ( 3+  sqrt( 3) )/  (4*sqrt(2));  beta4=(-l-sqrt(3))/  (4*sqrt(2)); 

%C  reate  1st  level  Wavelets,  then  from  this  generate  the  next  levels  in  the  M  RA 
W  ( 1, 1:4, 1)  =  [beta  I,beta2,beta3,beta4]; 
for  m=  2:( N /  2) 

W  (m,:,l)=  [W  (1,N  -2*(m-l)  +  1:N  ,1),W  (l,l:N-2*(m-l),l)  ]; 
end 

%C reate  2nd  through  10th  level  Waveletsfrom  1st  level  Wavelets 
for  j  =  2:10 

for  m  =  l:N/(2/vj)-l 

W  (m,:,j)=  betal*V (2*m-l,:,j-l)  +  beta2*V (2*m,:,j-l)  +  ... 
beta3*V  (2*m+  l,:,j-l)  +  beta4*V  (2*m+2,:,j-l); 
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end 

%Account  for  wrap  around  [W_  2m  +  l  &  W_  2m  +  2  will  exceed  N/(2~j)] 
m=N/(2/'j); 

W(m,:,j)  =  betal*V  (2*m-l,:,j-l)  +  beta2*V(2*m,:,j-l)  +  ... 

beta3*V  (2*m  +  l-N/ ( 2"(j -1) ) , :  ,j  - 1)  +  beta4*  V  ( 2*  m  +  2- N  /  ( 2"(j  - 1) ) ,:  ,j -1) ; 


Mathematica®  notebooks 


HSl.nb 

%T  his  notebook  ..nds  the  optic  fow  for  a  very  simple  approximation  to  motion, 

%  i.e.  motion  from  one  "image"  to  another  is  represented  as  a  bump  being 
%  added  to  the  surface  describing  the  grayscale  intensity  of  the  ..rst  image. 
<<Graphics'P  lotField'; 

%  T  he  parameter  nn  determines  how  many  modes  will  be 
%  used  in  the  approximation,  herenn  =  2. 
nn  =  2; 

%  T  his  creates  the  basis  for  approximating  u  and  v 
phi[i_  ,x_  ]=Cos[i  P  i  x]; 
u  =  Sum[a[i,j]phi[i,x]phi[j,y],{i,0,nn},{j,0,nn}]; 
v  =  Sum[b[i,j]phi[i,x]phi[j,y],{i,0,nn},{j  ,0,nn}]; 

%  T  his  creates  the  two  "images"  in  the  sequence,  el  and  e2, 

%  the  second  being  the  ..rst  with  a  bump  added 
%  to  represent  motion  of  grayscale  intensity. 
al=  0.45;bl=0.53; 

r[x_  ,y_  ,al_  ,bl_  ]=  Sqrt[(x-al) "2+ (y-bl) ^2]; 

el=Sin[3  Pi  y]Cos[2  Pi  x]; 

e2=el  +  lf[r[x,y ]<  1/  8,Cos[4  Pi  r[x,y, 0.45.0.53]], 0]; 

%  opl  and  op2  are  the  Euler-Lagrange  equations  where  the  parameter  eps  =  5. 
opl=D[el,x]/v2  u  +  D  [el,x]D[el,y]  v  -  eps  (D[u,{x,2}]+D  [u ,{y, 2}])  +  D[el,x](e2-el); 
op2=D[el,y]/v2  v  +  D[el,x]D[el,y]  u  -  eps  (D[v,{x,2}]+D[v,{y,2}])  +  D[el,y](e2-el); 
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%  T  his  sets  up  the  integral  equations  to  be  solved  by  the  G  alerkin  method 

Do[term  =  Expand[opl  phi [i ,x]ph i[j ,y ]];equati on l[i ,j]=M  ap[l  ntegrate[#  ,  {x,0,l},{y,0,l}]&,term]=  =  0; 

Print[Date[]];Print["i  =  ",  i,  ",  j  =  ",j],{i ,0,nn },{j ,0,nn}]; 

Do[term  =  Expand[op2  phi [i ,x]ph i[j ,y ]];equati on 2[i ,j ]=  M  ap[l  ntegrate[#  ,  {x,0,l},{y,0,l}]&,term]=  =  0; 

Print[Date[]];Print["i  =  ",  i,  ",  j  =  " ,j ],{i ,0, nn }, {j ,0, nn}]; 
coen  =  F  latten[T  ab  le[{a[i  ,j  ],  b[i  ,j  ]},{i  ,0,  nn},  {j  ,0,n  n}]]; 

eqns=F latten [J  oin |T able[equationl[i,j],{i,0,nn},{j,0,nn}],T ab le[equati on 2[i ,j ], {i , 0, nn }, {j ,0, nn }]]]; 

%  T  his  solves  the  integral  equations  and  plots  the  optic  tow  vector  ..eld 
eps=0.01; 

sol  =  Solve[eqns,coen  ]; 
uu[x_  ,y_  ]  =  First[u/.sol]; 
vv[x_  ,y_  ]  =  F  irst[v/.sol]; 

PlotVectorF  i  ei  d  [{uu  [x,y  ],vv[x,y  ]},  {x,0,l},{y,0, 1}] 

HSGalerkin.nb 

%  T  his  notebook  ..nds  the  optic  tow  for  the  R  ubik's  C  ube  sequence 
%  using  the  G al er kin  method  in  the  Calculus  of  Variations  approach 
<<Graphics'P  lotField' 

<<Graphics'Arrow' 

g=  I  mport  ["  C  :\\W  indows\\Desktop\\FI  or nSchunck\\rubikseq.gif"  ]; 

%  T  he  images  are  ..rst  interpolated  and  then  the  partial  derivatives 

%  are  taken  and  then  multiplied  to  form  the  terms  in  the  integral  equations 

el= L  i  stl  nterpolati  on[g[[l,  1,1]],  { (0, 1},  {0, 1}}]; 

e3= L  i  stl  nterpolati  on[g[[3, 1,1]],  { (0, 1},  {0, 1}}]; 

et[x_  ,y_  ]=e3[x,y]-el[x,y]; 

ex[x_  ,y_  ]=D[el[x,y],x]; 

ey[x_  ,y_  ]= D  [el[x,y ],y ]; 

ex2[x_  ,y_  ]=ex[x,y]~2; 

ey2[x_  ,y_  ]=  ey [x,y]^2; 

exey[x_  ,y_  ]= ex[x  ,y ]ey [x  ,y ]; 

exet[x_  ,y_  ]=ex[x,y]et[x,y]; 

eyet[x_  ,y_  ]=  ey [x,y ]et[x,y ]; 
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<<ex2.sav; 

<<ey2.sav; 

<<exdt.sav; 

<<eyet.sav; 

<<exey.sav; 

(*  E  X  2=  Tabl e[ex2[x,y  ], {x,0, 1, 1/  240} ,{y,0, 1,1/  256}]; 

E  Y  2=T able[N  [ey2[x,y]],{x, 0,1,1/  240}, {y, 0,1,1/  256}]; 

EX  E Y  =T able[N  [exey [x,y ]], {x,0, 1, 1/ 240}, {y,0, 1,1/  256}]; 

EX  ET  =Table[N  [exet[x,y]],{x, 0,1,1/  240}, {y, 0,1,1/  256}]; 

EY  ET  =Table[N  [eyet[x,y]],{x, 0,1,1/  240}, {y, 0,1, 1/256}];*) 

X  =T ab le[x, {x, 0,1, N  [1/  240]}];Y  =  Tabl e[y, {y,0, 1, N  [1/  256]}]; 

(*P  H I  [i_  ,j_  ]=T able[C os[i  Pi  x]Cos[|  Pi  y],{x, 0,1, N  [1/240]}, {y, 0,1, N  [1/256]}];*) 

P H I [i _  ,j_  ]=  Outer |T  imes,Cos[i  Pi  X],Cos[j  Pi  Y  ]]; 

%  T  he  cosine  basis  functions  are  set  up  to  use  17  modes  of  frequency  to  approximate 

%  the  optic  t  ow  vectors 

nn  =  16; ph i[i_  ,x_  ]=Cos[i  Pi  x]; 

u=Sum[a[i,j]phi[i,x]phi[j,y],{i,0,nn},{j,0,nn}]; 

v=Sum[b[i,j]phi[i,x]phi[j,y],{i,0,nn},{j,0,nn}]; 

%  T  he  cosine  approximations  for  the  terms  in  the  integral  equations  are  found 
(*Print[Date[]]; 

Do[coenexet[k,l]=l/240  1/  256  Plus  @@  F  latten [E X  ET 
P  H I  [k,  I  ]],  {k,0,  nn},  {1 , 0,  nn  }]; 

Print[Date[]]; 

Do[coeneyet[k,l]=  1/  240  1/  256  Plus  @@  F  latten[EY  ET 
P  H I  [k,l  ]],  {k,0,nn},  {1,0,  nn  }]; 

Print[Date[]]; 

Do[coenexey[k,l]=  1/240  1/  256  Plus  @@  F  latten  [EX  EY 
P  H I  [k,  I  ]],  {k,0,  nn},  {1 , 0,  nn  }]; 

Print[Date[]]; 

Do[coenex2[k,l]=  1/  240  1/  256  Plus  @@  Flatten[EX2 
P  H I  [k,  I  ]],  {k,0,  nn},  {1 , 0,  nn  }]; 


Print[Date[]]; 

Do[coeney2[k,l]=  1/  240  1/  256  Plus  @@  Flatten[EY2 
PHI  [k,l]],{k,0,nn},{l,0,nn}]; 

COEFFEXET=Table[coeHexdt[k,l],{kI0,nn},{lI0,nn}]; 

Save["COEF  F  EXET.sav",COEF  FEXET  ]; 

COEF  FEY  ET  =T ab le[coen eyet[k, I], {k, 0,n n },{ 1 ,0, nn }]; 

Save["COEF  F  EY  ET.sav",COEF  FEY  ET  ]; 

COEF  FEX  EY  =  Tabl e[coenexey[k, I ],{k,0, nn},{l,0,nn}]; 

Save["COEFFEXEY.sav"  ,COEFFEXEY]; 

COEF  FEX  2=T able[coenex2[k,l],{k,0,nn},{l,0,nn}]; 

Save["COEF  F  EX2.sav",COEF  FEX 2]; 

COEF  FEY  2=T able[coeney2[k,l],{k,0,nn},{l,0,nn}]; 

Save["COE F  F  EY  2.sav",COE F  FEY  2];*) 

«COEFFEX2.sav; 

«COEFFEY2.sav; 

«CO  E  F F  EX  E Y  .sav; 

«COEFFEXET.sav; 

<<CO  EFF  EY  ET  .sav; 

%  T  he  integral  for  a  multiplication  of  3  basis  functions  is  carried  out  in  advance 

%  so  that  it  does  not  have  to  be  computed  each  time  in  the  future 

tri  pi  e=T  ablefl  ntegrat  e[phi  [i,x]phi  [j  ,x]phi[k,x],{x,0,l}],{i,0,nn},{j,0,nn},{k,0,nn}]; 


%  T  he  integral  equations  are  set  up 

Do[equationl[k,l]=COEFFEXET[[k+l,l  +  l]]+Sum[COEFFEX2[[i  +  l,j  +  l]]a[m,n]triple[[i  + 
l,m+  l,k+l]]triple[[j  +  l,n+l,l+  l]],{i,0,nn},{j,0,nn},{m,0,nn},{n,0,nn}]+  Sum[ 

COEF  FEX  EY  [[i+  l,j +  l]]b[m,n]triple[[i+  I,m+l,k+l]]tripie[[j  +  l,n  +  l,l  +  l]],{i,0,nn}, 
{j,0,nn},{m,0,nn},{n,0,nn}]  + 

eps/4  ((k  P i ) ^2+ ( I  Pi)~2)  a[k,l]=  =  0,{k,0,nn},{l,0,nn}]; 

Print[Date[],",  equationl  is  done"]; 

Do[equation2[k,l]=COE  F  FEY  ET  [[k+  l,l  +  l]]+Sum[CO  EFF  EX  EY  [[i+  l,j  +  l]]a[m,n]triple[[ 
i  +  l,m+  l,k+l]]triple[[j  +  l,n+l,l+  l]],{i,0,nn},{j,0,nn},{m,0,nn},{n,0,nn}]+Sum[ 

COEF  FEY  2[[i  +  l,j  +  l]]b[m,n]triple[[i  +  l,m+l,k+  l]]triple[[j  +  l,n+  l,l  +  l]],{i,0,nn},{ 
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j,0,nn},{m,0,nn},{n,0,nn}]  + 

eps/4((k  P i)^2+  (I  Pi)~2)  b[k,l]=  =0,{k,0, nn }, {l,0,n n}]; 

P ri nt[D ate[], "  ,  equation2  is  done"]; 

coen  =  F  lattenfT  able[{a[i  ,j],b[i  ,j]},{i,0,nn},{j,0,nn}]]; 

eqns=F latten [J  oin |T able[equationl[i,j],{i,0,nn},{j,0,nn}],T ab lefequati on 2[i ,j ], {i , 0, nn }, {j , 0, nn }]]]; 

%  T  he  optic  fow  vectors  are  solved  for  using  S  =  .01 
eps=0.01; 

sol  =  Solve[eqns,coen  ]; 

Print[Date[],"  sol  is  done"]; 
uu[x_  ,y_  ]=  F i rstfu/ .sol ]; 
vv[x_  ,y_  ]=F irst[v/.sol]; 

discu=Table[uu[x,y],{x, 0,1,1/  240}, {y, 0,1,1/  256}]; 
discv=Table[vv[x,  y],  {x, 0,1,1/  240},  {y,  0,1, 1/  256}]; 

%  T  he  optic  fow  ..eld  is  plotted 

scale=  M  ax|T ab le[Sqrt[d i scu[[i ,j ]]^2  +  d iscv[[i ,j ]]^2], {i ,2,239},{j ,2, 255}]]; 
graph=  Show[T  able[G  raphics[A  rrow[{j,i  },{j  ,i  }+ 

100/ seal e{di scv[[i ,j ]],di scu[[i,j ]  ]},H eadScal in g- >R el ati ve],  Frame->T ruq 
P I ot L abel -> Stri ngj  oi n [" G al erkin  Method,  eps  =  ",  T oStringfeps],  ",  n  =  ",ToString[nn]]], 

{i, 2, 220, 20}, {j, 2, 256, 20}]] 


